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Abstract 



DDSCAT 7.2 is a freely available open-source Fortran-90 software package applying the "discrete 
dipole approximation" (DDA) to calculate scattering and absorption of electromagnetic waves by targets 
with arbitrary geometries and complex refractive index. The targets may be isolated entities (e.g., dust 
particles), but may also be 1-d or 2-d periodic arrays of "target unit cells", which can be used to study 
absorption, scattering, and electric fields around arrays of nanostructures. 

The DDA approximates the target by an array of polarizable points. The theory of the DDA and its 



implementation in DDSCAT is presented in Draine ( 1988 1 and Draine & Flatau ( 1994 1, and its extension 



to periodic structures in Draine & Flatau (2008 1. Efficient near-field calculations are carried out as 
described in |Flatau & Draine, ( |2012[ l. DDSCAT 7.2 allows accurate calculations of electromagnetic 
scattering from targets with "size parameters" 27raoff/A ^ 25 provided the refractive index m is not 
large compared to unity {\m — 1| ;S 2). DDSCAT 7.2 includes support for MP I, OpenMP, and the 
Intel® Math Kernel Library (MKL). 

DDSCAT supports calculations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, 
rectangular solids, finite cylinders, hexagonal prisms, etc.). Target materials may be both inhomoge- 
neous and anisotropic. It is straightforward for the user to "import" arbitrary target geometries into the 
code. DDSCAT automatically calculates total cross sections for absorption and scattering and selected 
elements of the Mueller scattering intensity matrix for specified orientation of the target relative to the 
incident wave, and for specified scattering directions. DDSCAT 7.2 can calculate scattering and absorp- 
tion by targets that are periodic in one or two dimensions. DDSCAT 7.2 calculates and stores the electric 
field E throughout a user-specified rectangular volume containing the target. A Fortran-90 code readnf 
to read E and P from files created by DDSCAT 7.2 nearfield calculations is included in the distribution. 

This User Guide explains how to use DDSCAT 7.2 (release 7.2.0) to carry out electromagnetic 
scattering calculations. 
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1 Introduction 

DDSCAT is a software package to calculate scattering and absorption of electromagnetic waves by tar- 
gets with arbitrary geometries using the "discrete dipole approximation" (DDA). In this approximation 
the target is replaced by an array of point dipoles (or, more precisely, polarizable points); the electromag- 
netic scattering problem for an incident periodic wave interacting with this array of point dipoles is then 
solved essentially exactly. The DDA (sometimes referred to as the "coupled dipole approximation") was 
apparently first proposed by Purcell & Pennypacker|( [T973] l. DDA theory was reviewed and developed 



further by |Draine| ( |1988| l, [Draine & Go odman (1993]l, reviewed by|Draine & Flatau|(il994) and |Draine 



( 2000 1, and recently extended to periodic structures by |Draine & Flatau ( 2008| l. 



DDSCAT 7.2, the current release of DDSCAT, is an open-source Fortran 90 implementation of the 
DDA developed by the authors [^DDSCAT 7.2 calculates absorption and scattering by isolated targets, or 
targets that are periodic in one or two dimensions, using methods described by |Draine & Flatau] ( |2008| l. 

DDSCAT is intended to be a versatile tool, suitable for a wide variety of applications including 
studies of interstellar dust, atmospheric aerosols, blood cells, marine microorganisms, and nanostructure 
arrays. As provided, DDSCAT 7.2 should be usable for many applications without modification, but the 
program is written in a modular form, so that modifications, if required, should be fairly straightforward. 

The authors make this code openly available to others, in the hope that it will prove a useful tool. We 
ask only that: 

• If you publish results obtained using DDSCAT, please acknowledge the source of the code, and 
cite relevant papers, such as IDraine & FlatauJ ( |1994| l, Draine & Flatau ( 2008 1, Flatau & Draine 
^2012) , and this UserGuide (Draine & Flatau 2012|l. 

• If you discover any errors in the code or documentation, please promptly communicate them to the 
authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public License) of 
the Free Software Foundation: you may copy, distribute, and/or modify the software identified as 
coming under this agreement. If you distribute copies of this software, you must give the recipients 
all the rights which you have. See the file doc/copylef t . txt distributed with the DDSCAT 
software. 

We also strongly encourage you to send email to draine@astro.princeton.edu identifying 
yourself as a user of DDSCAT; this will enable the authors to notify you of any bugs, corrections, or 
improvements in DDSCAT. Up-to-date information on DDSCAT can be found at 

| http : / / code . google . com/p/ddscat/| 
The latest version of DDSCAT 7.2 can be found at this URL. 

The current version, DDSCAT 7.2, uses the DDA formulae fromDraine|( 1988|, w ith dipole polariz- 



abihties determined from the Lattice Dispersion Relation ( Draine & Goodm an|1993[ Gutkowicz-Krusin 



|& Draine|2004| l. The code incorporates Fast Fourier Transform (FFT) methods (jCoodman etal.|1990F r 
DDSCAT 7.2 includes capability to calculate scattering and absorption by targets that are periodic in 



jusin] 
990r 



The release history of DDSCAT is as follows: 
DDSCAT 4b: Released 1993 March 12 
DDSCAT 4bl: Released 1993 July 9 

DDSCAT 4c: Although never announced, DDSCAT . 4 c was made available to a number of interested users beginning 1994 
December 18 

DDSCAT 5a7: Released 1996 
DDSCAT 5a8: Released 1997 April 24 
DDSCAT 5a9: Released 1998 December 15 
DDSCAT 5al0: Released 2000 June 15 
DDSCAT 6.0: Released 2003 September 2 
DDSCAT 6.1: Released 2004 September 10 
DDSCAT 7.0: Released 2008 September 1 
DDSCAT 7.1: Released 2010 February 7 
DDSCAT 7.2: Released 2012 February 15 



APPLICABILITY OF THE DDA 



one or two dimensions - arrays of nanostructures, for example. The theoretical basis for application of 
the DDA to periodic structures is developed in'Draine & Flatau (2008'). DDSCAT 7.2 includes capabil- 
ity to efficiently perform "nearfield" calculations of E in and around the target using FFT methods, as 



described by Flatau & Draine (2012 1. 

We refer you to the list of references at the end of this document for discussions of the theory and 
accuracy of the DDA [in particular, reviews by Draine & Flatau (19941 and Draine ( 2000] !, and recent 
extension to 1-d and 2-d arrays by [Draine & Flatau (2008 1.]. 

In ^we summarize the applicability of the DDA, and in ^we describe what the current release can 
calculate. 

In ^we describe the principal changes between DDSCAT 7.2 and the previous releases. The suc- 
ceeding sections contain instructions for: 

• obtaining the source code (^; 

• compiling and linking the code (^; 

• running a sample calculation (^; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations (Q; 

• specifying target orientation; 

• changing the DIMENSIONing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a Linux system; slight 
changes will be necessary for non-Linux sites, but they are quite minor and should present no difficulty. 

Finally, the current version of this User Guide can be obtained from 
fhttp: //arxiv.org/abs/12 02 . 342 4| 

Important Note: DDSCAT 7.2 differs in a number of respects from previous versions of DDSCAT. 
DDSCAT 7.2 includes support for both MPI and OpenMP, but - as of this writing - DDSCAT 7.2 has 
not yet been tested with either MPI or OpenMP. DDSCAT 7.2 has been tested extensively on single- 
processor systems, but if you are intending to use DDSCAT 7.2 with OpenMP or MPI, beware! Do at 
least a few comparison calculations in single-cpu mode to verify that the results with OpenMP or MPI 
appear to be correct. If you do encounter problems with OpenMP or MPI, please document them and 
communicate them to the authors. And if you find that everything appears to work properly, we'd like to 
know that too! 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, 
being limited only by the need to use an interdipole separation d small compared to (1) any structural 
lengths in the target, and (2) the wavelength A. Numerical studies ( [Draine & Goodman] 1993 Draine & 



|Flatau| 1 994[ |Draine|2000[ ) indicate that the second criterion is adequately satisfied if 

\m\kd < 1 , (1) 

where m is the complex refractive index of the target material, and /c = 27r/A, where A is the wavelength 
invacuo. This criterion is valid provide |to— 1| ;$ 3orso. When Im(TO) becomes large, the DDA solution 
tends to overestimate the absorption cross section Cabs, and it may be necessary to use interdipole 
separations d smaller than indicated by eq. ([T]) to reduce the errors in Cabs to acceptable values. 

If accurate calculations of the scattering phase function (e.g., radar or lidar cross sections) are desired, 
a more conservative criterion 

\m\kd < 0.5 (2) 

will usually ensure that differential scattering cross sections dCsca/d^ are accurate to within a few 
percent of the average differential scattering cross section Csca/47r ( Draine [2000 see). 



2 APPLICABILITY OF THE DDA 



8 



Let V be the actual volume of solid material in the target]^ If the target is represented by an array of 
N dipoles, located on a cubic lattice with lattice spacing d, then 



V = Ncf . 

We characterize the size of the target by the "effective radius" 



(3) 



(4) 



the radius of an equal volume sphere. A given scattering problem is then characterized by the dimen 
sionless "size parameter" 

2nacs 



A 



The size parameter can be related to N and \m\kd: 

27racff 62.04 / N 



A 



Equivalently, the target size can be written 



Oeff = 9.873- 



106 



1/3 



\m\kd 



A 



N 



ml V106 



1/3 



• \m\kd 



(5) 



(6) 



(7) 



Practical considerations of CPU speed and computer memory currently available on scientific worksta- 
tions typically limit the number of dipoles employed to < 10^ (see § 17 for limitations on N due to 
available RAM); for a given N, the limitations on |m|fc(i translate into limitations on the ratio of target 
size to wavelength. 

For calculations of total cross sections Cabs and Csca, we require \m\kd < 1: 



Ocff < 9.: 



A 



|m| V106 



N 



1/3 



or X < 



62.04 / N 



1/3 



For scattering phase function calculations, we require \m\kd < 0.5: 



aeft < 4.94t 



N 



1/3 



or a; < 



31.02 / N 



1/3 



(8) 



(9) 



It is therefore clear that the DDA is not suitable for very large values of the size parameter x, or very 
large values of the refractive index m. The primary utility of the DDA is for scattering by dielectric 
targets with sizes comparable to the wavelength. As discussed by Draine & Goodman ( 1993| l, Draine 
|& Flatau] ( |1994| l, and [Drairie] ( |2000| l, total cross sections calculated with the DDA are accurate to a few 
percent provided N > lO'' dipoles are used, criterion ([T]) is satisfied, and the refractive index is not too 
large. 

For fixed \m\kd, the accuracy of the approximation degrades with increasing |m — 1|, for reasons 



having to do with the surface polarization of the target as discussed by |ColUnge & Draine (20041. With 
the present code, good accuracy can be achieved for |to — 1| < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs.[T]-^ which show overall scattering 
and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a 
sphere, with N ranging from 304 to 59728. The DDA calculations assumed radiation incident along the 
(1,1,1) direction in the "target frame". Figs.[3}Qshow the scattering properties calculated with the DDA 
for X = ka — 7. Additional examples can be found in Draine & Flatau ( 1994| l and Draine (2000 1. 

As discussed below, DDSCAT 7.2 can also calculate scattering and absorption by targets that are 
periodic in one or two directions - for examples, see [Draine & Flatau] ( |2008| l. 



In the case of an infinite periodic target, V is the volue of soUd material in one "Target Unit Cell" 
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Figure 1: Scattering and absorption for a sphere with m = 1.33 + O.Oli. The upper panel shows the exact values of Qsca and 
Qabs, obtained with Mie theory, as functions of a; = ka. The middle and lower panels show fractional errors in Qsca and Qabs, 
obtained using DDSCAT with polarizabilities obtained from the Lattice Dispersion Relation, and labelled by the number A'^ of 
dipoles in each pseudosphere. After Fig. 1 of |Draine & Flatauljl994| l. 




3 DDSCAT 7.2 

3.1 What Does It Calculate? 

3.1.1 Absorption and Scattering by Finite Targets 

DDSCAT 7.2, like previous versions of DDSCAT, solves the problem of scattering and absorption by 
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on 

Figure 3: Differential scattering cross section for m — 1.33 + O.Oli pseudosphere and ka — 7. Lower panel shows fractional 
error compared to exact Mie theory result. The A'^ — 17904 pseudosphere has \m\kd — 0.57, and an rms fractional error in da/dQ 
of 2.4%. After Fig. 5 of |Draine & Flatau| ( [T994) . 




a finite target, represented by an array of polarizable point dipoles, interacting with a monochromatic 
plane wave incident from infinity. DDSCAT 7.2 has the capability of automatically generating dipole 
array representations for a variety of target geometries (see ^20] ) and can also accept dipole array repre- 
sentations of targets supplied by the user (although the dipoles must be located on a cubic lattice). The 
incident plane wave can have arbitrary elliptical polarization (see 5 23 1, and the target can be arbitrarily 
oriented relative to the incident radiation (see 5 18 i. The following quantities are calculated by DDSCAT 
7.2: 
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• Absorption efficiency factor Qabs = Cabs/Tra^ff, where Cabs is the absorption cross section; 

• Scattering efficiency factor Qsca = C'sca/'i'acfr' where Csca is the scattering cross section; 

• Extinction efficiency factor Qcxt = Qsca + Qahs', 

• Phase lag efficiency factor Qpha, defined so that the phase-lag (in radians) of a plane wave after 
propagating a distance L is just ntQphaTra^gi, where rit is the number density of targets. 

• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering properties of 



the target for scattering directions specified by the user (see 5 25 i. 
• Radiation force efficiency vector Qrad (see 5 16 1. 



• Radiation torque efficiency vector Qr (see ^16] ). 
In addition, the user can choose to have DDSCAT 7.2 store the solution for post-processing. 



3.1.2 Absorption and Scattering by Periodic Arrays of Finite Structures 

DDSCAT 7.2 includes the capability to solve the problem of scattering and absorption by an infinite 
target consisting of a 1-d or 2-d periodic array of finite structures, illuminated by an incident plane wave. 
The finite structures are themselves represented by arrays of point dipoles. 



The electromagnetic scattering problem is formulated by Draine & Flatau ( 2008 1, who show how the 
problem can be reduced to a finite system of linear equations, and solved by the same methods used for 
scattering by finite targets. 

The far-field scattering properties of the 1-d and 2-d periodic arrays can be conveniently represented 



by a generali zation of the Mueller scattering matrix to the 1-d or 2-d periodic geometry - see Draine & 



Flatau 



(|2008| for definition of S^^'^'' {M, () and Si-"' (M, N). For targets with 1-d periodicity, DDSCAT 



7.2 calculates S'|j''^(M, C) for user-specified M and For targets with 2-d periodicity, DDSCAT 7.2 

calculates (M, N) for both transmission and reflection, for user-specified (M, N). 

As for finite targets, the user can choose to have DDSCAT 7.2 store the calculated polarization field 
for post-processing. 



3.2 Application to Targets in Dielectric Media 

Let Lu be the angular frequency of the incident radiation. DDSCAT 7.2 now facilitates calculation of 
absorption and scattering by targets immersed in dielectric media (e.g., liquid water). In the parameter 
file ddscat . par, the user simply specifies the refractive index mmodium of the ambient medium. If the 
target is in vacuo, set mmcdium = 1- Otherwise, set mmedium to be the refractive index in the ambient 
medium at frequency lu. For example, H2O has mmedium — 1.335 near lu/2tt = 6 x 10^'* Hz, the 
frequency corresponding to Avac — 500 nm (green light). At this frequency, air at STP has mmedium = 
1.00028, which can be taken to be 1 for most applications. 

The wavelength provided in the file ddscat . par should be the vacuum wavelength Avac — 2ttc/uj 
corresponding to the frequency uj. 

The dielectric function or refractive index for the target material is provided via a file, with the 
filename provided via ddscat . par. The file should give either the actual complex dielectric function 
Ctargct or actual complex refractive index mtargct = Retarget of the target material, as a function of 
wavelength in vacuo. 

Internal to DDSCAT 7.2, the scattering calculation is carried out using the relative dielectric function 



relative refractive index: 



mrci(w) 



etargct('^) 
^medium ('^) 

mtargct(a^) 
^medium 



(10) 



(11) 
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and wavelength in the ambient medium 

^medium — ■ (12) 



''medium 



The absorption, scattering, extinction, and phase lag efficiency factors Qabs, Qsca, Qoxt, and Qpha 
calculated by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and 
extinction divided by Tra^g. For example, the attenuation coefficient for radiation propagating through a 
medium with a number density rit of scatterers will be just a = ntQcxtTra^^ . Similarly, the phase lag 
(in radians) after propagating a distance L will be utQpha'^algL. 

The elements 5^ of the 4x4 Mueller scattering matrix S calculated by DDSCAT for finite targets 
will be correct for scattering in the medium: 



A, 



medium 



2 



S • lin, (13) 



27rr 

where Ijn and Isca are the Stokes vectors for the incident and scattered light (in the medium), r is the 



distance from the target, and Amcdium is the wavelength in the medium (eq. 12 1. See 525 for a detailed 
discussion of the Mueller scattering matrix. 



The time-averaged radiative force and torque (see 5 16 1 on a finite target in a dielectric medium are 



Frad = QprTra^gUrad , (14) 
T-i r\ 2 Amcdium 

1 rad = t4r7raeffUrad , (15) 

ZTT 

where the time-averaged energy density is 

^lad f^mcdium o 7 v^^/ 

where Eq cos(ajt + 0) is the electric field of the incident plane wave in the medium. 

4 What's New? 

DDSCAT 7.2 differs from DDSCAT 7.1 in several ways, including: 

1. N.B.: The structure of the parameter file ddscat . par has again been changed: ddscat . par 
files that were used with DDSCAT 7.1 will need to be modified. See ^and Appendix [A] 

2. DDSCAT 7.2 now requires the user to specify the (real) refractive index rriambient of the ambient 
medium. See ^ 



DDSCAT 7.2 includes support for two additional conjugate gradient solvers (see 5 12 1: 

• GPBICG - Implementation by Chamuet and Rahmani of the Complex Conjugate Gradient 
solver presented by |Tang et al.| ( |2004) l. 

• QMRCCG - Quasi-Minimum-Residual Complex Conjugate Gradient solver, adapted from 



fortran-90 implementation kindly made available by RC. Chaumet and A. Rahmani ( ,Chaumet 
|& Rahmani|2009| l. 

DDSCAT 7.2 now allows the user to specify the maximum number of iterations (this was previ- 
ously hard-wired). 

DDSCAT 7.2 now supports fast calculations of the electric field within and near the target using 
FFT methods, as described by Flatau & Draine ( 2012] l; see ^28| ). The program DDFIELD that 



came with DDSCAT 7.1 is no longer needed (and no longer supported). 
The DDSCAT 7.2 distribution includes a set of "example" calculations: 
• FROM_FILE: target whose geometry is supplied via an ascii file 
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• ELLIPSOID: sphere, with nearfield calculation of E in and around the sphere. 

• RCTGLPRSM: rectangular brick 

• SPH_ANI_N: cluster of spheres, each characterized by an anisotropic dielectric function 

• CYLNDRPBC: infinite cylinder, calculated as a periodic array of disks. 

• DSKRCTPBC: doubly-periodic array of Au disks supported by a Si3N4 slab. 

• RCTGL_PBC: doubly-periodic array of rectangular blocks. 

• RCTGL_PBC: doubly-periodic array of rectangular blocks, with nearfield calculation of E 
in and around the sphere. 

• SPHRN_PBC: doubly-periodic array of clusters of spheres. 

• ELLIPSOID_NEARFIELD: same calculation as in the ELLIPSOID example, followed by 
a nearfield calculation of E in and around the sphere, followed by evaluation of E along a 
track passing through the sphere, plus creation of "VTK" files for visualization. 

• RCTGLPRSM_NEARFIELD: same calculation as in the RCTGLPRSM example, fol- 
lowed by a nearfield calculation of E in and around the target, followed by evaluation of 
E along a track passing through the prism, plus creation of "VTK" files for visualization. 

• RCTGL_PBC_NEARFIELD: same calculation as in the RCTGL_PBC example, followed 
by a nearfield calculation of E in and around the target, followed by evaluation of E along a 
track passing through the prism, plus creation of "VTK" files for visualization. 

7. DDSCAT 7.2 is distributed with a program VTRCONVERT.f90 that supports visualization of 
target geometries using the Visualization Toolkit (VTK), an open-source, freely-available software 
system for 3D computer graphics (http://www.vtk.org ) and, specifically. Para View ( ht tp://paraview.oi^ . 

8. DDSCAT 7.2 is distributed with a program READNF.f90 that 

• Allows the user to easily extract E at points along a line. 

• Creates "VTK" files for visuaUzation of |E| or |Ep using the VTK tools. 



5 Downloading the Source Code and Example Calculations 

DDSCAT 7.2 is written in standard Fortran-90 with a single extension: it uses the Fortran-95 standard 
library call CPU_TIME to obtain timing information. DDSCAT 7.2 is therefore portable to any system 
having a f90 or f95 compiler. It has been successfully compiled with many different compilers, including 
gfortran, g95, if ort, pgf 77, and NAG®f 95. 

It is possible to use DDSCAT 7.2 on PCs running Microsoft Windows operating systems, including 
Vista and Windows 7. Section[7]provides instructions for creating a unix-like environment in which you 
can compile and run DDSCAT 7.2. 



Alternatively, you may be able to obtain a precompiled native executable - see 5 7.1 More informa- 
tion on how to do this can be found at 

Ihttp : / /code . google . com/p/ ddscat| 

The remainder of this section will assume that the installation is taking place on a Unix, Linux or Mac 
OSX system with the standard developer tools (e.g., tar, make, and a f90 or f95 compiler) installed. 

The complete source code for DDSCAT 7.2 is provided in a single gzipped tarfile. To obtain the 
source code, simply point your browser to http : //code . google . com/p/ddscat/ and down- 
load the latest release. 

After downloading ddscat7 .2.0. tgz into the directory where you would like DDSCAT 7.2 to 
reside (you should have at least 10 Mbytes of disk space available), the source code can be installed as 
follows: 

If you are on a Linux system, you should be able to type 
tar xvzf ddscat7 . 2 . . tgz 
which wilF'extract" the files from the gzipped tarfile. If your version of "tar" doesn't support the "z" 
option, then try 

zcat ddscat7 . 2 . . tgz | tar xvf - 
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If neither of the above work on your system, try the two-stage procedure 

gunzip ddscat7 . 2 . . tgz 

tar xvf ddscatV . 2 . . tar 
The only disadvantage of the two-stage procedure is that it uses more disk space, since after the second 
step you will have the uncompressed tarfile ddscat? .2.0. tar - about 3.8 Mbytes - in addition to 
all the files you have extracted from the tarfile - another 4.6 Mbytes. 

Any of the above approaches should create subdirectories src, doc, diel, and examples_exp. 

• src contains the source code. 

• doc contains documentation (including this UserGuide). 

• diel contains a few sample files specifying refractive indices or dielectric functions as functions 
of vacuum wavelength. 

It is also reconmiended that you download ddscat? . 2 . 0_examples . tgz, followed by 
tar xvzf ddscat? . 2 . 0_examples . tgz 
Subdirectory examples_exp contains sample ddscat .par files as well as output files for various 
example problems, including both isolated targets and infinite periodic targets. It also includes sample 
readnf . par files for running readnf to support visuaUzation following nearfield calculations. 

6 Compiling and Linking on Unix/Linux Systems 

In the discussion below, it is assiuned that the DDSCAT 7.2 source code has been installed in a directory 
DDA/ src. The instructions below assume that you are on a Unix, Linux, or Mac OSX system. 

6.1 The default Makefile 

It is assumed that the system has the following already installed: 

• aFortran-90 compiler (e.g., gf ortran, g95, Intel® if ort, orNAG®f 95) . 

• cpp - the "C preprocessor". 

There are a number of different ways to create an executable, depending on what options the user wants: 

• what compiler and compiler flags? 

• single- or double-precision? 

• enable OpenMP? 

• enable MKL? 

• enable MPl? 

Each of the above choices needs requires setting of appropriate "flags" in the Makefile. 
The default Makefile has the following "vaniUa" settings: 

• gfortran -02 

• single-precision arithmetic 

• OpenMP not used 

• MKL not used 

• MPl not used. 

To compile the code with the default settings, simply position yourself in the directory DDA/ src, and 
type 

make ddscat 

If you have gfortran and cpp installed on your system, the above should work. You will get some 
warnings from the compiler, but the code will compile and create an executable ddscat. 

If you wish to use a different compiler (or compiler options) you will need to edit the file Makefile 
to change the choice of compiler (variable FC), compilation options (variable FFLAGS), and possibly 
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and loader options (variable LDFLAGS). The file Makefile as provided includes examples for com- 
pilers other than gf ortran; you may be able to simply "comment out" the section of Makefile that 
was designed for gf ortran, and "uncomment" a section designed for another compiler (e.g., Intel® 
if ort). 

6.2 Optimization 

The performance of DDSCAT 7.2 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags. 

6.3 Single vs. Double Precision 

DDSCAT 7.2 is written to allow the user to easily generate either single- or double -precision versions 
of the executable. For most purposes, the single-precision precision version of DDSCAT 7.2 should 
work fine, but if you encounter a scattering problem where the single-precision version of DDSCAT 7.2 
seems to be having trouble converging to a solution, you may wish to try using the double-precision 
version - this can be beneficial in the event that round-off error is compromising the performance of 
the conjugate-gradient solver. Of course, the double precision version will demand about twice as much 
memory as the single-precision version, and will take somewhat longer per iteration. 
The only change required is in the Makefile: for single-precision, set 

PRECISION = sp 
or for double-precision, set 

PRECISION = dp 

After changing the PRECISION variable in the Makefile (either sp dp, or dp sp), it is necessary 
to recompile the entire code. Simply type 

make clean 

make ddscat 
to create ddscat with the appropriate precision. 

6.4 OpenMP 

OpenMP is a standard for support of shared-memory parallel programming, and can provide a perfor- 
mance advantage when using DDSCAT 7.2 on platforms with multiple cpus or multiple cores. OpenMP 
is supported by many common compilers, such as gf ortran and Intel® if ort. 

If you are using a multi-cpu (or multi-core) system with OpenMP (www . openmp . org) installed, 
you can compile DDSCAT 7.2 with OpenMP directives to parallelize some of the calculations. To do 
so, simply change 

DOMP = 

OPENMP = 

to 

DOMP = -Dopenmp 
OPENMP = -openmp 
Note: OPENMP is compiler-dependent: gf ortran, for instance, requires 

OPENMP = -fopenmp 



6.5 MKL: the Intel® Math Kernel Library 

Intel® offers the Math Kernel Library (MKL) with the ifort compiler. This library includes DFTI for 
computing FFTs. At least on some systems, DFTI offers better performance than the GPFA package. 
To use the MKL library routine DFTI: 

• You must obtain the routine mkl_dfti . f 90 and place a copy in the directory where you are 
compiUng DDSCAT 7.2. mkl_dfti . f 90 is Intel® proprietary software, so we cannot distribute 
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it with the DDSCAT 7.2 source code, but it should be available on any system with a license for 
the Intel® MKL library. If you cannot find it, ask your system administrator. 

• Edit the Makefile: define variables CXFFTMKL . f and CXFFTMKL . o to: 

CXFFTMFKL.f = $(MKL_f) 
CXFFTMKL. f = ${MKL_o) 

• Successful linking will require that the appropriate MKL libraries be available, and that the string 
LFLAGS in the Makefile be defined so as to include these libraries. 

-Imk:l_em64t -lmkl_intel_thread -lmk:l_core -Iguide -Ipthread -lmk;l_intel_lp64 
appears to work on at least one installation. The Makefile contains examples. 

• type 

make clean 
make ddscat 

• The parameter file ddscat . par should have FFTMKL as the value of CMETHD. 



6.6 MPI: Message Passing Interface 

DDSCAT 7.2 includes support for parallelization under MPI. MPI (Message Passing Interface) is a 
standard for communication between processes. More than one implementation of MPI exists (e.g., 
mpich and openmpi). MPI support within DDSCAT 7.2 is compliant with the MPl-1.2 and MPI- 
2 standard^ and should be usable under any implementation of MPI that is compatible with those 
standards. 

Many scattering calculations will require multiple orientations of the target relative to the incident 
radiation. For DDSCAT 7.2, such calculations are "embarassingly parallel", because they are carried out 
essentially independently. DDSCAT 7.2 uses MPI so that scattering calculations at a single wavelength 
but for multiple orientations can be carried out in parallel, with the information for different orientations 
gathered together for averaging etc. by the master process. 

If you intend to use DDSCAT 7.2 for only a single orientation, MPI offers no advantage for DDSCAT 
7.2 so you should compile with MPI disabled. However, if you intend to carry out calculations for 
multiple orientations, and would like to do so in parallel over more than one cpu, and you have MPI 
installed on your platform, then you will want to compile DDSCAT 7.2 with MPI enabled. 

To compile with MPI disabled: in the Makefile, set 
DMPI = 

MIP . f = mpi_fake.f90 
MPI.o = mpi_fake.o 
To compile with MPI enabled: in the Makefile, set 

DMPI = -Dmpi 
MIP . f = $ (MPI_f ) 
MPI . o = $ (MPI_o) 

and edit LFLAGS as needed to make sure that you link to the appropriate MPI library (if in doubt, consult 
your systems administrator). The Makefile in the distribution includes some examples, but library 
names and locations are often system-dependent. Please do not direct questions regarding LFLAGS to 
the authors - ask your sys-admin or other experts familiar with your installation. 

Note that the MPI-capable executable can also be used for ordinary serial calculations using a single 
cpu. 



6.7 Device Numbers IDVOUT and IDVERR 

So far as we know, there are only one operating-system-dependent aspect of DDSCAT 7.2: the device 
number to use for "standard output". 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error mes- 
sages", respectively. Normally these would both be set to the device number for "standard output" (e.g., 

http : / /www . mpi-f orum. org/ 
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writing to the screen if running interactively). Variables IDVERR are set by DATA statements in the 
"main" program DDSCAT . f 90 and in the output routine WRIMSG (file wrimsg. f 90). The executable 
statement IDVOUT=0 initializes IDVOUT to 0. In the as-distributed version of DDSCAT . f 90, the state- 
ment 

OPEN (UNIT=IDVOUT, FILE=CFLLOG) 
causes the output to UNIT=IDVOUT to be buffered and written to the file ddscat . log_nnn, where 
nnn=0 for the first cpu, 001 for the second cpu, etc. If it is desirable to have this output unbuffered 
for debugging purposes, (so that the output will contain up-to-the-moment information) simply comment 
out this OPEN statement. 

7 Information for Windows Users 

There are several options to run DDSCAT on Windows. One can purchase a Fortran compiler such as 
Portland Group PGF compiler (http://www.pgroup.com/|. However, this option is not free. Below we 
discuss four methods which give access to DDSCAT on Windows using Open Source applications. We 
hav e tested each of these four methods. 



Method 


Compiler 


Advantage 


Problems 


Comments 


Native 


gfortran 


pre-compiled 


limiting options 


simplest 


MINGW 


gfortran 


compile with user options 


compilation step 


simple 


Cygwin 


gfortran, G95 


compile with user options 


compilation step 


simple 


Virtualbox/UBUNTU 


gfortran, G95, Intel 


full LINUX access 


learning curve 


difficult 



7.1 Native executable 

We provide a self-extracting executable which includes a pre-compiled, ready-to-run, DDSCAT exe- 
cutable. The advantage is that user avoids the possibly difficult compilation step, and has immediate ac- 
cess to DDSCAT. Beginning with release 7.2 we package the Windows distribution using "Inno Setup5" 
(|http://www.jrsoftware.org ) which is a free installer for Windows programs. This will install a Windows 
native executalbe ddscat .exe as well as source code, documentation, and relevant test examples. 
This is by far the the simplest way to get DDSCAT running on a Windows system. However, we provide 
only single precision version without optimization and with limitation of target size. 

Thus, serious DDSCAT users may, at some stage, need to recompile the code as outlined below. 
The executable should be executed using the windows "cmd" command which opens a separate shell 
window. One then has to change directory to where DDSCAT is placed, and execute DDSCAT invoking 

ddscat.exe 

7.2 Compilation using MINGW 

The executable file discussed in "native executable" section was compiled using "gfortran" in the MINGW 
environment (http://www.mingw.org/). Windows users may choose to recompile code using MINGW. 
MINGW and MSYS provide several tools crucial for compilation of DDSCAT to a native windows 
executable. MinGW ("Minimalist GNU for Windows") is a minimalist development environment for na- 
tive Microsoft Windows applications. It provides a complete Open Source programming tool set which 
is suitable for the development of native MS-Windows applications, and which do not depend on any 
3rd-party libraries. 

MSYS ("Minimal SYStem"), is a Bourne Shell command line interpreter system. Offered as an 
alternative to Microsoft's cmd.exe, this provides a general purpose command line environment, which 
is particularly suited to use with MinGW, for porting applications to the MS-Windows platform. We 
suggest using the installer mingw-get-inst (for example mingw-get-inst-20111118 . exe) 
which is a simple Graphical User Interface installer that installs MinGW and MSYS. During the GUI 
phase of the installation select the default options, but from the following list you need to select Fortran 
Compiler, MSYS Basic System and MinGW Developer ToolKit: 

1 . MinGW Compiler Suite C Compiler optional 
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2. C++ Compiler 

3. optional Fortran Compiler 

4. optional ObjC Compiler 

5. optional MSYS Basic System 

6. optional MinGW Developer Toolkit 

You will have to add PATH to bin directories as described in the MINGW. 

Once you have MINGW and MSYS installed, in "all programs" there is now the MinGW program 
"MinGW shell". Once you open this shell one can now compile Fortran code using "gfortran". The 
"make" and "tar" utilities are available. The command "pwd" shows which directory corresponds to the 
initial "MinGW" shell. For example it can be "/home/Piotr" which corresponds to windows directory 
"c:\mingw\msys\l .O\home\piotr". 

You now need to copy ddscat. tar files and untar them. To make a native executable, edit "Makefile" 
so that the "LFLAGS" string is defined to be 

LFLAGS=-static-libgcc -static -libgfortran 

and then execute "make all". The resulting ddscat . exe doesn't require any non-windows libraries. 
This can be checked with the 

objdump -X ddscat.exe I grep "DLL" 
command. 

7.3 Compilation using CYGWIN 

Another option is provided by "CYGWIN", an easily installable UNIX-like emulation package. It is 
available from http://www.cygwin.com/ . It installs authomatically. However, during installation one has 
to specify installation of several packages incuding the f95 compiler gfortran, the make utility, the tar 
utility, and nano. Once installed, you will be able to open the CYGWIN shell and make DDSCAT using 
the standard Linux commands as discussed in this manual (see 

7.4 UBUNTU and Virtualbox 

By far the most comprehensive solution to running DDSCAT on windows is to install UBUNTU Linux 
under Oracle Virtualbox. First install Oracle Virtualbox (from https://www.virtualbox.org/) and then 
install UBUNTU Linux (from http://www.ubuntu.com/). You will be able to run a full LINUX environ- 
ment on your Windows computer You can add (using "synaptic file manager") gfortran and many other 
packages including graphics. Once a f90 or f95 compiler has been installed, you can compile DDSCAT 
as described in ^ 

8 A Sample Calculation: RCTGLPRSM 

When the tarfile is unpacked, it will create four directories: src, doc, diel, and examples_exp. 
The examples_exp directory has a number of subdirectories, each with files for a sample calculation. 
Here we focus on the files in the subdirectory RCTGLPRSM. To follow this, go to the DDA directory (the 
directory in which you unpacked the tarfile) and 

cd examples_exp 
to enter the examples_exp directory, and 

cd RCTGLPRSM 
to enter the RCTGLPRSM directory. 

9 The Parameter File ddscat . par 

It is assumed that you are positioned in the examples_exp/RCTGLPRSM. directory, as per ^ The 
file ddscat . par (see also Appendix [A]) provides parameters to the program ddscat: 
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' Parameter file for v7 . 2 ===================' 

'■k-k-k* Preliminaries 

'NOTORQ' = CMT0RQ*6 (NOTORQ, DOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, GPBICG, PETRKP, QMRCCG) — solution method 
'GPFAFT' = CMDFFT*6 (GPFAFT, FFTMKL) — FFT method 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) — prescription for polari zabi 1 it ies 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) — specify binary output 
'■k*** Initial Memory Allocation ****' 

100 100 100 = dimensioning allowance for target generation 
' Target Geometry and Composition kkkk' 
'RCTGLPRSM' = CSHAPE*9 shape directive 
16 32 32 = shape parameters 1-3 

1 = NCOMP = number of dielectric materials 
' . . /diel/Au_evap' - file with refractive index 1 
'k*** Additional Nearfield calculation? 

= NRFLD (=0 to skip nearfield calc, =1 to calculate nearfield E) 

0.0 0.0 0.0 0.0 0.0 0.0 (tract, extens . of calc. vol. in -x, +x, -y , +y , -z , +z ) 
'kkkk Error Tolerance kkkk' 

l.OOe-5 - TOL - MAX ALLOWED (NORM OF | G>-AC| E>-ACA | X>) / (NORM OF AC|E>) 

maximum number of iterations allowed kkkk' 
300 - MXITER 

'kkkk Interaction cutoff parameter for PBC calculations kkkk' 
l.OOe-2 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

Angular resolution for calculation of <cos>, etc. kkkk' 
0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] ^^2 ) 

Vacuum wavelengths (micron) kkkk' 
0.5000 0.5000 1 'LIN' = wavelengths ( fir st , last , how many , how=LIN, INV, LOG) 
'kkkk Refractive index of ambient medium' 
1.000 - NAMBIENT 

Effective Radii (micron) kkkk ' 
0.246186 0.246186 1 'LIN' = aeff ( first , last , how many , how=LIN, INV, LOG) 
'kkkk Define Incident Polarizations kkkk' 

(0,0) (l.,0.) (0.,0.) - Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol. state eOl; =2 to also do orth. pol. state) 
'kkkk Specify which output files to write kkkk' 

1 = IWRKSC (-0 to suppress, -1 to write ".sea" file for each target orient. 

Prescribe Target Rotations kkkk' 
0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 0. 1 - THETMI, THETMX, NTHETA (theta-angle between al and k) 

0. 0. 1 - PHIMIN, PHIMAX, NPHI (phi-rotation angle of al around k) 

'kkkk Specify first IWAV, IRAD, lORI (normally 0) kkkk' 
= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'kkkk Select Elements of S_ij Matrix to Print kkkk' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions kkkk' 

'LFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

2 = NPLANES = number of scattering planes 

0. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 1 
90. 0. 180. 5 - phi, thetan_min, thetan_max, dtheta (in deg) for plane 2 

Here we discuss the general structure of ddscat . par (see also Appendix [A|). 



9.1 Preliminaries 

ddscat . par starts by setting the values of five strings: 

• CMDTRQ specifying whether or not radiative torques are to be calculated (e.g., NOTORQ) 

• CMDSOL specifying the solution method (e.g., PBCGS2) (see section JT2]i. 



• CMDFFT specifying the FFT method (e.g., GPFAFT) (see section ^ i. 

• CALPHA specifying the polarizability prescription (e.g., GKDLDR) (see § 14 1. 

• CBINFLAG specifying whether to write out binary files (e.g., NOTBIN) 
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9.2 Initial Memory Allocation 

Three integers 

MXNX MXNY MXNZ 

are given that need to be equal to or larger than the anticipated target size (in lattice units) in the xtf, 
Ytf, and ztf directions. This is required only for initial memory allocation for the target-generation 
stage of the calculation. For example, 
100 100 100 

would require initial memory allocation of only ~ 100 MBytes. Note that after the target geometry has 
been determined, DDSCAT 7.2 wiU proceed to reallocate as much memory as is actually required for 
the calculation. 

9.3 Target Geometry and Composition 

• CSHAPE specifies the target geometry (e.g., RCTGLPRSM) 

As provided, the file ddscat . par is set up to calculate scattering by a 32x64x64 rectangular array of 
131072 dipoles. 

The user must specify NCOMP, the number of different compositions (i.e., dielectric functions) that 
will be used. This is then followed by NCOMP Unes, with each line giving the name (in quotes) of a 
dielectric function file. In our example, NCOMP is set to 1. The dielectric function of the target material 
is provided in the file . . / diel/Au_evap, which gives the refractive index of evaporated Au over a 
range of wavelengths. The file Au_evap is located in subdirectory examples_exp/diel . 

9.4 Additional Nearfield Calculation? 

The user set NRFLD = or 1 to indicate whether the first DDSCAT calculation (solving for the target 
polarization, and absorption and scattering cross sections) should automatically be followed by a second 
"nearfield" calculation to evaluate the electric field E throughout a rectangular volume containing the 
original target. 

• If NRFLD = , the nearfield calculation will not be done. 
The next line in ddscat . par wiU be read but not used. 

• If NRFLD = 1 , the nearfield calculation will be done. 

The next line in ddscat . par then specifies 6 non-negative numbers, ri, r2, ra, r^, r^, re specifying the 
increase in size of the nearfield computational volume relative to the circumscribing rectangular volume 
used for the original solution. If the original volume is Xi < x < X2, Yi < y < Y2, Z\ < z < Z^, then 
the nearfield calculation will be done in a volume {Xi — tyL^) < x < {X2 + ^2^2.), (Yi — r-^Ly) < 
V < iY2 + TiLy), {Zi - rsL,) < z < {Z2 + reL,), where = {X2 - Xi), Ly = {Y2 - Yi), 
L, = {Z2 - Zi). 

9.5 Error Tolerance 

TOL = the error tolerance. Conjugate gradient iteration will proceed until the linear equations are solved 
to a fractional error TOL. The sample calculation has TOL=10~^. 

9.6 Maximum Number of Iterations 

MX ITER = maximum number of complex-conjugale-gradient iterations allowed. As there are 3A^ equa- 
tions to solve, MX ITER should never be larger than 3A^, but in practice should be much smaller. As 
a default we suggest setting MXITER=100, but for some problems that converge very slowly you may 
need to use a larger value of MX ITER. 
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9.7 Interaction Cutoff Parameter for PBC Calculations 

GAMMA = parameter limiting certain summations that are required for periodic targets ( [Draine & Flatau| 



2008 see). The value of GAMMA has no effect on computations for finite targets - it can be set to any 

value, including 0). 

For targets that are periodic in 1 or 2 dimensions, GAMMA needs to be small for high accuracy, but 
the required cpu time increases as GAMMA becomes smaller. GAMMA = 10^^ is reasonable for initial 
calculations, but you may want to experiment to see if the results you are interested in are sensitive to 
the value of GAMMA. pDraine & Flatau|p008| l show examples of how computed results for scattering can 
depend on the value of 7. 



9.8 Angular Resolution for Computation of (cos 9), etc. 

The parameter ETASCA determines the selection of scattering angles used for computation of certain 
angular averages, such as {cos 6) and (cos^0), and the radiation pressure force (see 5 16 1 and radia- 
tive torque (if CMDTRQ=DOTORQ). Small values of ETASCA result in increased accuracy but also cost 
additional computation time. ETASCA=0.5 generally gives accurate results. 

If accurate computation of (cos 9) or the radiation pressure force is not required, the user can set 
ETASCA to some large number, e.g. 10, to minimize unnecessary computation. 



9.9 Vacuum Wavelengths 

Wavelengths A (in vacuo) are specified in one line in ddscat . par consisting of values for 4 variables: 

WAVINI WAVEND NWAV CDIVID 
where WAVINI and WAVEND are real numbers, NWAV is an integer, and CDIVID is a character variable. 

• If CDIVID = ' LIN' , the A will be uniformly spaced beween WAVINI and WAVEND. 

• If CDIVID = ' INV , the A will be uniformly spaced in 1/A between WAVINI and WAVEND 

• If CDIVID = ' LOG' , the A will be uniformly spaced in log(A) between WAVINI and WAVEND 

• If CDIVID = ' TAB' , the A will be read from a user-supplied file wave . tab, with one wave- 
length per line. For this case, the values of WAVINI, WAVEND, and NWAV will be disregarded. 

The sample ddscat . par file specifies that the calculations be done for a single wavelength (A = 0.50). 
The units must be the same as the wavelength units used in the file specifying the refractive index. 
In this case, we are using /zm. 

9.10 Refractive Index of Ambient Medium, mmedium 

In some cases the target of interest will be immersed in a transparent medium (e.g., water) with a re- 
fractive index different from vacuum. The refractive index TOmcdium should be specified here. As the 
medium is assumed to be transparent, mmedium is a real number. DDSCAT will calculate the scatter- 
ing properties for the target immersed in the ambient medium. If the target is located in a vaccum, set 

mnicdium 1- 



9.11 Target Size aes 

Note that in DDSCAT the "effective radius" acfr is the radius of a sphere of equal volume - i.e., a sphere 
of volume Nd'^ , where d is the lattice spacing and N is the number of occupied (i.e., non-vacuum) lattice 
sites in the target. Thus the effective radius Ccff = (3Af/47r)^/'^(i . Our target should have a thickness 
a = 0.5fj,m in the xtf direction. If the rectangular solid is a x b x c, with a : 6 : c :: 32 : 64 : 64, then 
V = abc = 4a3. Thus a^s = {W/Att)^/^ = {3/Tr)^^^a = 0.49237/im. 

The target sizes Ocff to be studied are specified on one line of ddscat . par consisting of 4 variables: 
AEFFINI AEFFEND NRAD CDIVID 
where tt AEFFINI and AEFFEND are real numbers, NRAD is an integer, and CDIVID is a character 
variable. 
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• If CDIVID = ' LIN' , the flcff will be uniformly spaced beween AEFFINI and AEFFEND. 

• IfCDIVID = ' INV, the flcff will be uniformly spaced in 1/acff between AEFFINI and AEFFEND 

• If CDIVID = ' LOG', the floff will be uniformly spaced in log(aoff) between AEFFINI and 
AEFFEND 

• IfCDIVID = ' TAB' , the Ocff will be read from a user-supplied file aeff . tab, with one value 
of Oeff per line, beginning on line 2. For this case, the values of AEFFINI, AEFFEND, and NRAD 
will be disregarded (DDSCAT will read and use all the values of floff in the table (beginning on 
line 2). 

The sample ddscat .par file specifies that the calculations be done for a single acff = 0.49237/im. 
The target is a rectangular solid with aspect ratio 1:2:2. The thickness in the x direction is (7r/3)^/'^acff = 
0.500/im. 



9.12 Incident Polarization 

The incident radiation is always assumed to propagate along the xlf axis - the a;-axis in the "Lab 
Frame". The sample ddscat . par file specifies incident polarization state Gqi to be along the ylf axis 
(and consequently polarization state eo2 will automatically be taken to be along the zlf axis). I0RTH=2 
in ddscat .par calls for calculations to be carried out for both incident polarization states (gqi and 
eo2 - see ^23). 



9.13 Target Orientation 



The target is assumed to have two vectors ai and a2 embedded in it; a2 is perpendicular to ai . For the 
present target shape RCTGLPRSM, the vector ai is along the xtf axis of the target, and the vector a2 is 



along the yxF axis (see §20.18 i. The target orientation in the Lab Frame is set by three angles: j3, 9, and 
defined and discussed below in 5 18 Briefly, the polar angles 6 and $ specify the direction of ai in 



the Lab Frame. The target is assumed to be rotated around ai by an angle /?. The sample ddscat . par 
file specifies /? = and $ = (see lines in ddscat .par specifying variables BETA and PHI), and 
calls for three values of the angle 8 (see line in ddscat . par specifying variable THETA). DDSCAT 
chooses Q values uniformly spaced in cos 9. In this case we specify 1 : we obtain only one value: 
9 = 0.0 



9.14 Starting Values of I WAV, IRAD, lORl 

Normally we begin the calculation with IWAVE=0, IRAD=0, and IORI=0. However, under some cir- 
cumstances, a prior calculation may have completed some of the cases. If so, the user can specify starting 
values of IWAV, IRAD, TORI; the computations will begin with this case, and continue. 



9.15 Which Mueller Matrix Elements? 

The sample parameter file specifies that DDSCAT should calculate 6 distinct Mueller scattering matrix 
elements Sij, with 11, 12, 21, 22, 31, 41 being the chosen values of ij. 



9.16 What Scattering Directions? 
9.16.1 Isolated Finite Targets 

For finite targets, the user may specify the scattering directions in either the Lab Frame (' LFRAME' ) or 
the Target Frame (' T FRAME' ). 

For finite targets, such as specified in this sample ddscat .par, the Sij are to be calculated for 
scattering directions specified by angles (6*, 0). 

''Had we specified 90 3 we would have obtained three values of O between and 90°, unformly-spaced in cos 6: Q = 0, 60°, 
and 90°. 
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The sample ddscat .par specifies that 2 scattering planes are to be specified: the first has (j) = 
and the second has (p = 90°; for each scattering plane 9 values run from to 180° in increments of 10°. 



9.16.2 1-D Periodic Targets 

For periodic targets, the scattering directions must be specified in the Target Frame: ' TFRAME' 

CMDFRM . 

Scattering from 1-d periodic targets is discussed in detail by 'Draine & Flatau ( 2008| l. For periodic 



targets, the user does not specify scattering planes. For 1-dimensional targets, the user specifies scattering 
cones, corresponding to different scattering orders M. For each scattering cone M, the user specifies 
Cmin, Cmax, (in degrees); the azimuthal angle ( will run from (^in to (max, in increments of A^. For 
example: 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 
1 = number of scattering cones 

0. 0. 180. 0.05 = OrderM zetamin zetamax dzeta for scattering cone 1 



9.16.3 2-D Periodic Targets 

For targets that are periodic in 2 dimensions, the scattering directions must be specified in the Target 
Frame: ' TFRAME' = CMDFRM. 



Scattering from targets that are periodic in 2 dimensions is discussed in detail by Draine & Flatau 



(2008 L For 2-D periodic targets, the user specifies the diffraction orders (Af, N) for transmitted radi- 



ation: the code will automatically calculate the scattering matrix elements S'^J'^'' (M, N) for both trans- 
mitted and reflected radiation for each (Af , N) specified by the user. For example: 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 

0. 0. = OrderM OrderN for scattered radiation 



10 Running DDSCAT 7.2 Using the Sample ddscat. par File 

It is again assumed that you are in directory . . /DDA/examples_exp/RCTGLPRSM, as per ^ The 
ddscat executable (created as per the instructions in ^ is assumed to be . . /DDA/ src/ddscat 



10.1 Single-Process Execution 

To execute the program on a UNIX system (running either sh or csh), simply create a symbolic link by 
typing 

In -s ../../src/ddscat ddscat 
or you could simply move the previously-created executable into the current directory (assumed to be 
../DDA/examples_exp/RCTGLPRSM/) by typing 

mv ../../src/ddscat ddscat 
Then, to perform the calculation, type 

ddscat >& ddscat. out & 
which will redirect the "standard output" to the file ddscat . out, and run the calculation in the back- 
ground. 

The sample calculation [32x64x64=131072 dipole target, 3 target orientations, two incident polar- 
izations for each orientation, with scattering (Mueller matrix elements Sij) calculated for 37 distinct 
scattering directions], requires 672 cpu sec to complete on a 2.53 GHz cpu. for the assumed Au compo- 
sition, between 28 and 32 iterations were required for the complex conjugate gradient solver to converge 
to the specified error tolerance of TOL = 1 . e-5 for each orientation and incident polarization. 
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10.2 Code Execution Under MPI 

Local installations of MPI will vary - you should consult with someone familiar with way MPI is in- 
stalled and used on your system. 

At Princeton University Dept. of Astrophysical Sciences we use PBS (Portable Batch System^ to 
schedule jobs. MPI jobs are submitted using PBS by first creating a shell script such as the following 
example file pbs . submit: 

#!^in/bash 

#PBS -1 nodes=2:ppn=l 

#PBS -1 mem=1200MB,pmem=300MB 

#PBS -m bea 

#PBS -j oe 

cd $PBS_0_WORKDIR 
/usr/local/bin/mpiexec ddscat 

The lines beginning with #PBS -1 specify the required resources: 

#PBS -1 nodes=2 : ppn=l specifies that 2 nodes are to be used, with 1 processor per node. 
#PBS -1 mem=1200MB, pmem= 3 MB specifies that the total memory required (mem) is 1 200MB, 
and the maximum physical memory used by any single process (pmem) is 300MB. The actual definition 
of mem is not clear, but in practice it seems that it should be set equal to 2x(nodes)x(ppn)x(pmem) . 
#PBS -m bea specifies that PBS should send email when the job begins (b), and when it ends (e) or 
aborts (a). 

#PBS - j oe specifies that the output from stdout and stderr will be merged, intermixed, as stdout. 
This example assumes that the executable dscat is located in the same directory where the code is to 
execute and write its output. If ddscat is located in another directory, simply give the full pathname, 
to it. The qsub command is used to submit the PBS job: 

qsub pbs. submit 

As the calculation proceeds, the usual output files will be written to this directory: for each wave- 
length, target size, and target orientation, there will be a file viaaarbbb\ccc .sea, where aaa=000, 001, 
002, ... specifies the wavelength, bbb=000, 001, 002, ... specifies the target size, and cc=000, 001, 002, ... 
specifies the orientation. For each wavelength and target size there will also be a file waaarbbbori . avg 
with orientationally-averaged quantities. Finally, there will also be tables qtable, and qtable2 with 
orientationally-averaged cross sections for each wavelength and target size. 

In addition, each processor employed will write to its own log file ddscat . log_nnn, where 
nnn=000, 001, 002, .... These files contain information concerning cpu time consumed by different 
parts of the calculation, convergence to the specified error tolerance, etc. If you are uncertain about how 
the calculation proceeded, examination of these log files is recommended. 

11 Output Files 
11.1 ASCII files 

If you run DDSCAT using the command 
ddscat >& ddscat.log & 
you will have various types of ASCII files when the computation is complete: 

• afileddscat.log; 

• afilemtable; 

• a file qtable; 

' http : / /www . openpbs . org 
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• a file qtable2; 

• files vixxxryyyori. . avg (one, wOOOrOOOori. avg, for the sample calculation); 

• if ddscat . par specified IWRKSC=1, there will also be files vxxxryyy'Kzzz . sea (1 for the sam- 
ple calculation: wOOOrOOOkOOO . scalcl, wOOOrOOOkOOl . sea, w000r000k002 . sea). 

The file ddseat . out will contain minimal information (it may in fact be empty). 

The file ddseat . log_000 will contain any error messages generated as well as a running report 
on the progress of the calculation, including creation of the target dipole array. During the iterative 
calculations, Qoxt^ Qabs. and Qpha are printed after each iteration; you will be able to judge the degree 
to which convergence has been achieved. Unless TIME IT has been disabled, there will also be timing 
information. If the MP I option is used to run the code on multiple cpus, there will be one file of the form 
ddseat . log_nnn for each of the cpus, with nnn=000 , 001,002, . . .. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Qcxt> Qabs> Qsca> 
5(1) = (cos(6's)), (cos^(6's)), Qbk, and iVsca- Here Qoxt, Qabs, and Qsca are the extinction, absorption, 
and scattering cross sections divided by Tra^g. Qbk is the differential cross section for backscattering 
(area per sr) divided by tto^jj. A'sca is the number of scattering directions used for averaging over 
scattering directions (to obtain (cos 6*), etc.) (see 524i. 

The file qtable2 contains a summary of the orientationally-averaged values of Qpha. Qpoh and 
Qcpoi- Here Qpha is the "phase shift" cross section divided by ira^g ( Draine|[l988 see definition in). 
Qpoi is the "polarization efficiency factor", equal to the difference between Qext for the two orthogonal 
polarization states. We define a "circular polarization efficiency factor" Qcpoi = QpoiQpha, since an 
optically-thin medium with a small twist in the alignment direction will produce circular polarization in 
initially unpolarized light in proportion to Qcpoi- 

For each wavelength and size, DDSCAT 7.2 produces a file with a name of the form 
wxxjcryyyori . avg, where index xxx (=000, 001, 002....) designates the wavelength and index yyy 
(=000, 001, 002...) designates the "radius"; this file contains Q values and scattering information aver- 
aged over however many target orientations have been specified (see ^TS] The file wOOOrOOOori.avg 
produced by the sample calculation is provided below in Appendix [B] 

In addition, if ddseat .par has specified IWRKSC=1 (as for the sample calculation), DDSCAT 
7.2 will generate files with names of the form vxxxryyy'kzzz . avg, where xxx and yyy are as before, and 
index zzz =(000,001,002...) designates the target orientation; these files contain Q values and scattering 
information for each of the target orientations. The structure of each of these files is very similar to 
that of the vi xxxr yyy ori. . avg files. Because these files may not be of particular interest, and take up 
disk space, you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the 
sample calculation with IWRKSC=1. 

The sample ddseat . par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, 
and averaging over 3 target orientations. Running DDSCAT 7.2 with the sample ddseat . par file will 
therefore generate files wOOOrOOOkOOO . sea, wOOOrOOOkOOl . sea, and w000r000k002 . sea . 
To understand the information contained in one of these files, please consult Appendix[C] which contains 
an example of the file wOOOrOOOkOOO. sea produced in the sample calculation. 



11.2 Binary Option 

It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete information, 
including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in 
ddseat . par . 

Subroutine writebin . f 90 provides an example of how this can be done. The "header" section 
contains dimensioning and other variables which do not change with wavelength, particle geometry, 
and target orientation. The header section contains data defining the particle shape, wavelengths, parti- 
cle sizes, and target orientations. If ALLBIN has been specified, the "data" section contains, for each 
orientation, Mueller matrix results for each scattering direction. The data output is limited to actual di- 
mensions of arrays; e.g. nscat ,4,4 elements of Muller matrix are written rather than mxseat ,4,4. 
This is an important consideration when writing postprocessing codes. 
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12 Choice of Iterative Algorithm 

As discussed elsewhere (e.g., |Draine|[l988[ |Draine & Flatau||l994| l, the problem of electromagnetic 
scattering of an incident wave by an array of N point dipoles can be cast in the form 

AP = E (17) 

where E is a 3A^-dimensional (complex) vector of the incident electric field Einc at the N lattice sites, 
P is a 3A^-dimensional (complex) vector of the (unknown) dipole polarizations, and A is a 3N x 3A^ 
complex matrix. 

Because 3iV is a large number, direct methods for solving this system of equations for the unknown 
vector P are impractical, but iterative methods are useful: we begin with a guess (typically, P = 0) for 



the unknown polarization vector, and then iteratively improve the estimate for P until equation ( 17 1 is 
solved to some error criterion. The error tolerance may be specified as 

lAtAP - AtEl , 

^ <h , (18) 



|AtE 

where A^^ is the Hermitian conjugate of A [{A^)ij = (Aji)*], and h is the error tolerance. We typically 



use h — 10^^ in order to satisfy eq.( 17 1 to high accuracy. The error tolerance h can be specified by the 
user through the parameter TOL in the parameter file ddscat . par (see Appendix [A|. 

A major change in going from DDSCAT.4b to 5a (and subsequent versions) was the implementation 
of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT 
7.2 is now structured to permit solution algorithms to be treated in a fairly "modular" fashion, facilitating 
the testing of different algorithms. A number of algorithms were compared by Flatau (1997 p| two of 
them (PBCGST and PETRKP) performed well and were made available to the user in the DDSCAT 
6.0 release. DDSCAT 7.1 introduced a third option, PBCGS2. DDSCAT 7.2 includes two more CCG 
options: GPBICG and QMRCCG. The choice of algorithm is made by specifying one of the options (here 
in alphabetical order): 



GPBICG - Generalized Product-type methods based on Bi-CG Zhang] ( [T997 1. We use an imple 



mentation suggested by Tang et al. ( 2004)), and coded by P.C. Chaumet and A. Rahmani (Chaumet 



[&~R ahmani 2009]). We are grateful to P.C. Chaumet and A. Rahmani for making their code avail- 
able. 

PBCGS2 - BiConjugate Gradient with Stabilization as implemented in the routine ZBCG2 by 
M.A. Botchev, University of Twente. This is based on the PhD thesis of D.R. Fokkema, and on 
work by |Sleijpen & van der VOTSt] ( |1 995 1 [19961 1. 



PBCGST - Preconditioned BiConjugate Gradient with STabilitization method from the Parallel 
Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

PETRKP - the complex conjugate gradient algorithm of |Petravic & Kuo-Petravic| ( |1979[ ), as coded 
in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algo- 
rithm discussed by Draine ( 1988) 1 and used in the earliest versions of DDSCAT. 



• QMRCCG - the quasi-minimum-residual complex conjugate gradient algorithm, based on f77 code 
written by P.C. Chaumet and A. Rahmani, converted here to f90 and adapted to single/double 
precision. 

All five methods work fairly well. Our experience suggests that PBCGS2 is generally fastest and best- 
behaved, and we recommend that the user try it first. There have been claims that QMRCCG and/or 
GPBICG are faster, but this has not been our experience. PETRKP is slow but may prove stable on 
some problems where more aggressive algorithms become unstable. We have not carried out systematic 
studies of the relative performance of the different algorithms - the user is encouraged to experiment. 

If a fast algorithm the case it runs into numerical difficulties, PBCGST and PETRKP are available as 
alternatives^] 

* A postscript copy of this report - file eg . ps ~ is distributed with the DDSCAT 7.2 documentation. 

' The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rdd@ukc . ac . uk) and Tim Hopkins (trhgukc . ac . uk) 
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13 Choice of FFT Algorithm 

DDSCAT 7.2 offe rs two FFT options: (1) the GPFA FFT algorithm developed by Dr. Clive Temperton 
(|Temperton| 1992|), p]and (2) the Intel® MKL routine DFTI. 




Figure 5: Comparison of cpu time required by 3 different FFT implementations. It is seen that the GPFA and FFTW implemen- 
tations have comparable speeds, much faster than Brenner's FFT implementation. 



The GPFA routine is portable and quite fast: Figure 5 compares the speed of three FFT implemen- 
tations: Brenner's, GPFA, and FFTW (http : / /www . f f tw . org). We see that while for some cases 
FFTW 2 . 1 . 5 is faster than the GPFA algorithm, the difference is only marginal. The FFTW code and 
GPFA code are quite comparable in performance - for some cases the GPFA code is faster, for other 
cases the FFTW code is faster. For target dimensions which are factorizable as 2*3^5''" (for integer i, 
j, k), the GPFA and FFTW codes have the same memory requirements. For targets with extents A^^;, 

is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and scalar computers using a variety 
of iterative methods (available at 

,http : / /chasqueweb . uf rgs .br/~rudnei . cunha/pim. htmlll. PIM offers a number of iterative methods, including 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB ( [Van der Vorst|1992| l, 

• the restarted version of BICGSTAB, RBICGSTAB jSleijpen & Fokkema| ( [T993) 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can be called direcdy via 
ddscat.par. It is possible to add other options by changing the code in getfml.f90 . |Flatau| (1997| | has compared 
the convergence rates of a number of different methods. A helpful introduction to conjugate gradient methods is provided 
by the report "Conjugate Gradient Method Without Agonizing Pain" by lonathan R. Shewchuk, available as a postscript file: 
£tp: //RE PORTS .ADM. CS ■CMU.EDU/usr0/anon/1994/CMU-CS-94-125 .ps^ 

°The GPFA code contains a parameter LVR which is set in data statements in the routines gpf a2f, gpf a3f , and gpf a5f . 
LVR is supposed to be optimized to correspond to the "length of a vector register" on vector machines. As delivered, this parameter 
is set to 64, which is supposed to be appropriate for Crays other than the C90. For the C90, 128 is supposed to be preferable (and 
perhaps "preferable" should be read as "necessary" - there is some basis for fearing that results computed on a C90 with LVR other 
than 128 run the risk of being incorrect!) The value of LVR is not critical for scalar machines, as long as it is fairly large. We found 
little difference between LVR=64 and 128 on a Sparc 10/51, on an Ultrasparc 170, and on an Intel® Xeon cpu. You may wish to 
experiment with different LVR values on your computer architecture. To change LVR, you need to edit gpf a . f 90 and change the 
three data statements where LVR is set. 
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Figure 6: Scattering and absorption for a m = 1.7 + O.li sphere, calculated using two prescriptions for the polarizability: 
LATTDR is the lattice dispersion relation result of jDraine & Goodman| ( [T993) . GKDLDR is the lattice dispersion relation result 
of Gutkowicz-Krusin & Draine 1 2004 1. Results are shown as a function of scattering parameter x — 27ra / A; the upper scale gives 



values of |m|fc(i. We see that the cross sections calculated with these two prescriptions are quite similar for \m\kd ^ 0.5. For other 
examples see |Gutkowicz-Krusin & Draine| ( [2004} . 



Ny, Nz which are not factorizable as 2'3"'5'^, the GPFA code needs to "extend" the computational vol- 
ume to have values of N^, Ny, and which are factorizable by 2, 3, and 5. For these cases, GPFA 
requires somewhat more memory than FFTW. However, the fractional difference in required memory is 
not large, since integers factorizable as 2*3^ 5*^ occur fairly frequently]^ [Note: This "extension" of the 
target volume occurs automatically and is transparent to the user] 

DDSCAT 7.2 offers a new FFT option: the Intel® Math Kernel Library DFTI. This is tuned for 
optimum performance, and appears to offer real performance advantages on modern multi-core cpus. 
With this now available, the FFTW option, which had been included in DDSCAT 6.1, has been removed 
from DDSCAT 7.2. 

The choice of FFT implementation is obtained by specifying one of: 

• FFTMKL to use the Intel® MKL routine DFTI (see ^6.5| l. This is recommended, but requires 
that the Intel® Math Kernel Library be installed on your system. 

• GPFAFT to use the GPFA algorithm ( |Temperton|1992| i. This is not quite as fast as FFTMKL, 
but is written in plain Fortran-90. It is a perfectly good alternative if the Intel® Math Kernel 
Library is not available on your system. 



14 Dipole Polarizabilities 



Option GKDLDR specifies that the polarizability be prescribed by the "Lattice Dispersion Relation" 



with the polarizability found by Gutkowicz-Krusin & Draine (^WA ), who corrected a subtle error in the 



^2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 
125, 128, 135, 144, 150, 160, 162, 180, 192, 200 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 
432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 
1080, 1125, 1152, 1200, 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 1920, 1944, 2000, 
2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 
3600, 3645, 3750, 3840, 3888, 4000, 4050, 4096 are the integers < 4096 which are of the form 2*3-' 5'= . 
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analysis of Draine & Goodman ( 1993 k For \m\kd ^ 1, the GKDLDR polarizability differs slightly from 
the LATTDR polarizability, but the differences in calculated scattering cross sections are relatively small, 
as can be seen from Figure[6] We recommend option GKDLDR. 

Users wishing to compare can invoke option LATTDR to specify that the "Lattice Dispersion Rela- 
tion" of 'Draine & Goodman ( 1993 1 be employed to determine the dipole polarizabilities. This polariz- 
abihty also works well. 



15 Dielectric Functions 

In order to assign the appropriate dipole polarizabilities, DDSCAT 7.2 must be given the refractive index 
m or dielectric constant e of the material (or materials) of which the target of interest is composed. This 
information is supplied to DDSCAT 7.2 through a table (or tables), read by subroutine DIELEC in file 
dielec . f 90, and providing either the complex refractive index m = n + ik or complex dielectric 
function e = ei + ie2 as a function of wavelength A. Since m = e^/^, or e = iri^, the user must supply 
either m or e. 

DDSCAT 7.2 can calculate scattering and absorption by targets with anisotropic dielectric functions, 
with arbitrary orientation of the optical axes relative to the target shape. See |27| 

The table containing the dielectric function information should give to or e as a function of the 
wavelength in vacuo. 

The table formatting is intended to be quite flexible. The first line of the table consists of text, up 
to 80 characters of which will be read and included in the output to identify the choice of dielectric 
function. (For the sample problem, it consists of simply the statement m = 1.33 + O.Oli.) The 
second line consists of 5 integers; either the second and third or the fourth and fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(m) is stored in. 

• The third integer specifies which column Im(TO) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 

• The fifth integer specifies which column Im(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and Im(e) from the file; if the 
fourth and fifth integers are zeros, then Re(TO) and Im(m) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There must be 
at least 3 lines of data: even if to or e is required at only one wavelength, please supply two additional 
"dummy" wavelength entries in the table so that the interpolation apparatus will not be confused. 



As discussed in 5 3.2 DDSCAT can scattering for targets embedded in dielectric media. The re- 
fractive index of the ambient medium is specified by the value of NAMBIENT in the parameter file 
ddscat . par (see §9.10| l. 

Here is an example of a refractive index file for Au: 

Gold, evaporated (Johnson & Christy 1972, PRB 6, 4370) 
12 3 0= columns for wave. Re (n) , Im(n), epsl, eps2 



wave (um) 


Re (n) 


Im (n) 


epsl 


eps2 


0.5486 





43 


2 


455 


-5 


84 


2.11 


0.5209 





62 


2 


081 


-3 


95 


2.58 


0.4959 


1 


04 


1 


833 


-2 


28 


3.81 


0.4714 


1 


31 


1 


849 


-1 


70 


4 .84 


0.4509 


1 


38 


1 


914 


-1 


76 


5.28 


0.4305 


1 


45 


1 


948 


-1 


69 


5.65 


.4133 


1 


46 


1 


958 


-1 


70 


5 .72 


. 3974 


1 


47 


1 


952 


-1 


65 


5.74 


0.3815 


1 


46 


1 


933 


-1 


60 


5.64 


0.3679 


1 


48 


1 


895 


-1 


40 


5.61 


0.3542 


1 


50 


1 


866 


-1 


23 


5.60 


0.3425 


1 


48 


1 


871 


-1 


31 


5.54 



CALCULATION OF (COS 6), RADIATIVE FORCE, AND RADIATION TORQUE 



0.3315 1. 

0.3204 1. 

0.3107 1. 

0.3009 1. 



48 1.883 

54 1.898 

53 1.893 

53 1.889 



-1.36 5.57 

-1.23 5.85 

-1.24 5.79 

-1.23 5.78 



16 Calculation of (cos 6), Radiative Force, and Radiation Torque 

In addition to solving the scattering problem for a dipole array, DDSCAT can compute the three- 
dimensional force Fiad and torque Fiad exerted on this array by the incident and scattered radiation 
fields. The radiation torque calculation is carried out, after solving the scattering problem, only if 
DOTORQ has been specified in ddscat .par. For each incident polarization mode, the results are 
given in terms of dimensionless efficiency vectors Qpi and Qr, defined by 

Qp, = ^^ , (19) 

Qr^^J^^ , (20) 

where Fj-ad and F^ad are the time-averaged force and torque on the dipole array, k ~ 27r/A is the 
wavenumber in vacuo, and Uiad = Eq/Stt is the time-averaged energy density for an incident plane 
wave with amplitude Eq cos{ujt + 0). The radiation pressure efficiency vector can be written 

Qpr = Qcxtk - Qscag , (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean direction 
of propagation of the scattered radiation: 

1 f ^^ dC (n,k) ^ ^22^ 



^sca J 

where dil is the element of solid angle in scattering direction n, and dCsca/dO, is the differential scat- 
tering cross section. The components of Qpi are reported in the Target Frame: Qpr.i = Fjad • xtf, 

Qpr, 2 = Frad ' yXF, Qpr, 3 = Frad ' ^TF- 



Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner 
( |199 6). It is important to note that evaluation of Qpr and Qr involves averaging over scattering directions 
to evaluate the linear and angular momentum transport by the scattered wave. This evaluation requires 



appropriate choices of the parameter ETASCA - see §24 

In addition, DDSCAT calculates (cos 6) [the first component of the vector g in eq. ( 22 1] and the 
second moment (cos^ 0) . These two moments are useful measures of the anisotropy of the scattering. 
For example, |Draine| ( |2003[ ) gives an analytic approximation to the scattering phase function of dust 
mixtures that is parameterized by the two moments (cos 9) and (cos^ 0) . 



17 Memory Requirements 

The memory requirements are determined by the size of the "computational volume" - this is a rectan- 
gular region, of size NXxNYxNZ that is large enough to contain the target. If using the GPFAFT option, 
then NX, NY, NZ are also required to have only 2,3, and 5 as prime factors (see footnote|9|. 
In single precision, the memory requirement for DDSCAT 7.2 is approximately 

(35. + 0.0010 X NX X NY X NZ) Mbytes for single precision (23) 

(42 + 0.0020 X NX X NY X NZ) Mbytes for double precision (24) 
Thus, in single precision, a 48x48x48 calculation requires ^--^146 MBytes. 
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The memory is allocated dynamically - once the target has been created, DDSCAT 7.2 will deter- 
mine just how much overall memory is needed, and will allocate it. However, the user must provide 
information (via ddscat . par) to allow DDSCAT 7.2 to allocate sufficiently large arrays to carry out 
the initial target creation. Initially, the only arrays that will be allocated are those related to the target 
geometry, so it is OK to be quite generous in this initial allowance, as the memory required for the target 
generation step is small compared to the memory required to carry out the full scattering calculation. 



18 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x direction. 
For purposes of discussion we will always let unit vectors xlf, Ylf, zlf = xlf x Ylf be the three 
coordinate axes of the LF. 

In ddscat . par one specifies the first polarization state eoi (which obviously must lie in the y, z 
plane in the LF); DDSCAT automatically constructs a second polarization state 692 = xlf x orthog- 
onal to Gqi. Users will often find it convenient to let polarization vectors Gqi = y, eo2 = z (although 
this is not mandatory - see ^23| ). 




Fi gure 7 ! Target orientation in the Lab Frame, x = xlf is the direction of propagation of the incident radiation, and y = y lf 
is the direction of the real component (at xi^y = 0, t = 0) of the first incident polarization mode. In this coordinate system, the 
orientation of target axis ai is specified by angles and "I>. With target axis ai fixed, the orientation of target axis a2 is then 
determined by angle /9 specifying rotation of the target around ai. When /? = 0, a2 lies in the ai,XLF plane. 



Recall that definition of a target involves specifying two unit vectors, ai and a2, which are imagined 
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to be "frozen" into the target. We require a2 to be orthogonal to ai . Therefore we may define a "Target 
Frame" (TF) defined by the three unit vectors ai, 3.2, and 3.3 = ai x a2 . 

For example, when DDSCAT creates a 32x24x16 rectangular solid, it fixes ai to be along the 
longest dimension of the soUd, and a2 to be along the next-longest dimension. 

Important Note: for periodic targets, DDSCAT 7.2 requires that the periodic target have ai = xtf 
and a2 = ytf- 

Orientation of the target relative to the incident radiation can in principle be determined two ways: 

1 . specifying the direction of ai and a2 in the LF, or 

2. specifying the directions of xlf (incidence direction) and ylf in the TF. 

DDSCAT uses method 1.: the angles O, $, and /3 are specified in the file ddscat . par. The target is 
oriented such that the polar angles 8 and $ specify the direction of ai relative to the incident direction 
Xlf, where the xlf^Ylf plane has $ = 0. Once the direction of ai is specified, the angle /3 then 
specifies how the target is to rotated around the axis ai to fully specify its orientation. A more extended 
and precise explanation follows; 



18.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 9, and /3, to specify the directions of unit vectors ai and a2 in the LF 
(see Fig.|7]l. 

Q is the angle between ai and xlf- 

When $ = 0, ai will lie in the Xlf , Ylf plane. When $ is nonzero, it will refer to the rotation of ai 
around xlf: e.g., $ = 90° puts ai in the xlf; ^lf plane. 

When /? = 0, 3.2 will lie in the xlf, ai plane, in such a way that when 6 = and $ = 0, 3.2 is in 
the yLF direction: e.g, Q — 90°, ^ — Q, (3 — has 3i = ylf and 3.2 = — xlf- Nonzero f3 introduces 
an additional rotation of a2 around ai: e.g., Q ~ 90°, $ = 0, = 90° has ai ~ yLF and a2 = zlf- 

Mathematically: 

ai = Xlf cos 9 + ylf sin 9 cos $ + zlf sin 6 sin $ (25) 
812 — —Xlf sin 9 cos /3 + ylf [cos 9 cos (3 cos $ — sin /3 sin $] 

+Zlf [cos 9 cos /3 sin $ + sin (3 cos $] (26) 
3.3 — Xlf sin 9 sin f3 — jlf [cos 9 sin /3 cos $ + cos (3 sin $] 

— zlf[cos 9 sin/3 sin $ — cos /3 cos $] (27) 



or, equivalently: 



Xlf ~ ai cos 9 — 32 sin 9 cos (3 + sin 9 sin (3 (28) 
YLF = ai sin 9 cos $ + a2 [cos 9 cos /3 cos $ — sin (3 sin $] 

— a3[cos9sin/3cos$ + cos/3sin<i>] (29) 
Zlf = ai sin 9 sin $ + 3.2 [cos 9 cos f3 sin $ + sin /3 cos $] 

— 33[cos 9 sin /3 sin $ — cos /? cos <1>] (30) 



18.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that xlf (the direction 
of propagation of the radiation) and ylf (usually the first polarization direction) and zlf (= xlf x 
Ylf) refer to certain directions in the TF. Given the definitions of the LF and TF above, this is simply 
an exercise in coordinate transformation. For example, one might wish to have the incident radiation 
propagating along the (1,1,1) direction in the TF (example 14 below). Here we provide some selected 
examples: 

1. XLF = ai, YLF = a2, Zlf = a3 : 9 = 0, $ + /3 

2. XLF = ai, YLF = a3, zlf = -a2 : 9 = 0, $ + /3 = 90° 
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3. 


xlf 


= a2, Ylf = 


ai, Zlf = -as : 6 = 90°, P = 180°, $ = 


4. 


xlf 


= 32, YlF = 


as, Zlf = ai : e = 90°, P = 180°, $ = 90° 


5. 


Xlf 


= as, ylf = 


ai, Zlf = a2 : e = 90°, P = -90°, $ = 


6. 


Xlf 


= as, ylf = 


a2, Zlf = -ai : 6 = 90°, P = -90°, $ = -90° 


7. 


Xlf 


= -ai, Ylf 


= £12, Zlf = -as : 9 = 180°, P + <^ = 180° 


8. 


Xlf 


= -ai, ylf 


= as, Zlf = a2 : e = 180°, /3 + $ = 90° 


9. 


Xlf 


= -a2, Ylf 


= ai, Zlf = as : e = 90°, /3 = 0, = 


10. 


Xlf 


= -a2, Ylf 


= as, Zlf = -ai : 6 = 90°, /3 = 0, $ = -90° 


11. 


Xlf 


= -as, Ylf 


= ai, Zlf = -a2 : 6 = 90°, P = -90°, $ = 


12. 


Xlf 


= -as, Ylf 


= a2, Zlf = ai : e = 90°, P = -90°, $ = 90° 


13. 


Xlf 


= (ai + a2)/^/2, ylf = as, zlf = (ai - a2)/y2 : 6 = 45°, p 


14. 


Xlf 


= (ai + as - 


1- as)/y3, ylf = (ai - a2)/\/2, zlf = (ai + a2 - 




e = 


54.7356°, P 


= 135°, $ = 30°. 



90° 



18.3 Sampling in O, and /? 

The present version, DDSCAT7.2, chooses the angles P, 6, and $ to sample the intervals (BETAMI, BETAMX), 

(THETMI, THETMX) , (PHIMIN, PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, 
PHIMAX are specified in ddscat . par . The prescription for choosing the angles is to; 

• uniformly sample in P; 

• uniformly sample in 

• uniformly sample in cos 6. 

This prescription is appropriate for random orientation of the target, within the specified Umits of P, 
ande. 

Note that when DDSCAT 7.2 chooses angles it handles P and $ differently from 6. The range for p 
is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take BETAMI=0, 
BETAMX=90, NBETA=2 you will get ^ = 22.5° and 67.5°. Similarly, if you take PHIMIN=0, PHIMAX=180, 
NPHI=2 you will get $ = 45° and 135°. 

Sampling in Q is done quite differently from sampling in P and First, as already mentioned above, 
DDSCAT 7.2 samples uniformly in cos Q, not Q. Secondly, the samphng depends on whether NTHETA 
is even or odd. 

• If NTHETA is odd, then the values of 6 selected include the extreme values THETMI and THETMX; 
thus, THETMI=0, THETMX=90, NTHETA=3 will give you 6 = 0, 60°, 90°. 

• If NTHETA is even, then the range of cos 6 will be divided into NTHETA intervals, and the mid- 
point of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you 
e = 41.41° and 75.52° [cos 9 = 0.25 and 0.75]. 

The reason for this is that if odd NTHETA is specified, then the "integration" over cos 9 is performed 
using Simpson's rule for greater accuracy. If even NTHETA is specified, then the integration over cos 9 
is performed by simply taking the average of the results for the different 9 values. 

If averaging over orientations is desired, it is reconnmended that the user specify an odd value of 
NTHETA so that Simpson's rule will be employed. 



19 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about 
this depends on the distribution of orientations which is applicable. 
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19.1 Randomly-Oriented Targets 

For randomly-oriented targets, we wish to compute the orientational average of a quantity Q{P, 6, $): 

iQ) = ^l d$g(/3,e,$) . (31) 

To compute such averages, all you need to do is edit the file ddscat .par so that DDSCAT knows 
what ranges of the angles /S, O, and $ are of interest. For a randomly-oriented target with no symmetry, 
you would need to let /3 run from to 360°, 6 from to 180°, and $ from to 360°. 

For targets with symmetry, on the other hand, the ranges of (3, O, and $ may be reduced. First 
of all, remember that averaging over $ is relatively "inexpensive", so when in doubt average over to 
360°; most of the computational "cost" is associated with the number of different values of i(3,Q) which 
are used. Consider a cube, for example, with axis ai normal to one of the cube faces; for this cube /3 
need run only from to 90°, since the cube has fourfold symmetry for rotations around the axis ai. 
Furthermore, the angle 9 need run only from to 90°, since the orientation (/3,6,$) is indistinguishable 
from iP, 180° - 6, 360° - 

For targets with symmetry, the user is encouraged to test the significance of i3,Q,^ on targets with 
small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

19.2 Nonrandomly-Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around the x axis = 
direction of propagation of the incident radiation), one may be able to use DDSCAT 7.2 with appropriate 
choices of input parameters. More generally, however, you will need to modify subroutine ORIENT to 
generate a list of NBETA values of 0, NTHETA values of 6, and NPHI values of plus two weight- 
ing arrays WGTA ( 1-NTHETA, 1-NPHI ) and WGTB ( 1-NBETA) . Here WGTA gives the weights which 
should be attached to each (Q,^) orientation, and WGTB gives the weight to be attached to each /3 ori- 
entation. Thus each orientation of the target is to be weighted by the factor WGTAxWGTB. For the case 
of random orientations, DDSCAT chooses 9 values which are uniformly spaced in cos 9, and /3 and $ 
values which are uniformly spaced, and therefore uses uniform weights 
WGTB=1./NBETA 

When NTHETA is even, DDSCAT sets 
WGTA=1./(NTHETAXNPHI) 

but when NTHETA is odd, DDSCAT uses Simpson's rule when integrating over 9 and 
WGTA= (1/3 or 4/3 or 2/3)/(NTHETAxNPHI) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented 
cases. If, for example, the orientation is such that for a given $ value only one 9 value is possible 
(this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the 
Earth's atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the 
combinations of 9 and $ which the present version of DDSCAT is set up to do. We hope to improve 
this in a future version of DDSCAT. 

20 Target Generation: Isolated Finite Targets 

DDSCAT contains routines to generate dipole arrays representing finite targets of various geometries, 
including spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two touching 
ellipsoids, and three touching ellipsoids. The target type is specified by variable C SHAPE on line 9 of 
ddscat .par, up to 12 target shape parameters (SHPARi, SHPAR2, SHPAR3, ...) on line 10. The 
target geometry is most conveniently described in a coordinate system attached to the target which we 
refer to as the "Target Frame" (TF), with orthonormal unit vectors xtf, Ytf, ^tf = ^tf x Ytf- Once 
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the target is generated, the orientation of the target in the Lab Frame is accomphshed as described in 
T81 

Every target generation routine will specify 
The "occupied" lattice sites; 

The composition associated with each occupied lattice site; 

Two "target axes" ai and a2 that are used as references when specifying the target orientation; and 
The location of the Target Frame origin of coordinates. 
;et geometries currently supported include: 
FROM_FILE : isotropic target material(s), geometry read from file ( ^20.1 1 



Tari 



ANIFRMFIL : anisotropic target material(s), geometry read from file (i 20.2 1 



ANIELLIPS : anisotropic ellipsoid ( pOJ] ) 

ANI_ELL_2 : two touching anisotropic ellipsoids (single composition) (5 20.41 
ANI_ELL_3 : three touching anisotropic ellipsoids (single composition) ( ^20.5 i 
ANIRCTNGL : anisotropic brick ( pOg i 
CONELLIPS : two concentric ellipsoids) ( pOT i 
CYLINDERl : finite cylinder ( pOS i 

CYLNDRCAP : cyhnder with hemispherical end-caps (P0.9|) 



DSKRCTNGL : disk resting on a brick (5 20.10 1 



DW1996TAR : 13-block target used by Draine & Weingartner ( 1996 1 (5 20.11 



ELLIPSOID : ellipsoid (including spheroid and sphere) (5 20. 12 1 
ELLIPSO_2 : two touching ellipsoids, different compositions allowed (i 20.13 i 



ELLIPSO_3 : three touching ellipsoids, different compositions allowed (i 20. 14 1 



HEX_PRISM : finite hexagonal prism ( pTTS i 
LAYRDSLAB : multilayer rectangular slab ( ^20.16[ ) 
MLTBLOCKS : collection of cubic blocks (^2017} 



RCTGLPRSM : rectangular prism (i.e., brick) ( pOTS i 
RCTGLBLK3 : stack of 3 rectangular blocks ( p0.19[ ) 
SLAB_HOLE : rectangular slab with cylindrical hole ( |20.20 i 
SPHERES_N : collection of N spheres ( pOlT] ) 



SPHROID_2 : two touching spheroids, different compositions allowed (^ 20.22| i 
SPH_ANI_N : collection of N anisotropic spheres (520.23 i 
TETRAHDRN : tetrahedron ( |2a24] i 
TRNGLPRSM : triangular prism (^025 1 



UNIAXICYL : finite cylinder of uniaxial material ( ^20.26 1 
Each is described below. 



20.1 FROM_FILE = Target composed of possibly anisotropic material, defined 
by list of dipole locations and "compositions" obtained from a file 

If anisotropic, the "microcrystals" in the target are assumed to be aligned with the principal axes of 
the dielectric tensor parallel to xtf, Vtf, and ztf- This option causes DDSCAT to read the target 
geometry and composition information from a file shape . dat instead of automatically generating 
one of the geometries for which DDSCAT has built-in target generation capability. The shape . dat 
file is read by routine REASHP (file reashp.f90). The file shape.dat gives the number N of 
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dipoles in the target, the components of the "target axes" ai and a2 in the Target Frame (TF), the vector 
xo(l — 3) determining the correspondence between the integers IXYZ and actual coordinates in the TF, 
and specifications for the location and "composition" of each dipole. The user can customize REASHP as 
needed to conform to the manner in which the target description is stored in file shape . dat. However, 
as supplied, REASHP expects the file shape . dat to have the following structure: 

• one fine containing a description; the first 67 characters will be read and printed in various output 
statements 

• N = number of dipoles in target 

• dix o-iy o-iz = x,y,z components (in TF) of ai 

• a2x o,2y 0'2z = x,y,z Components (in TF) of a2 

• dx/d dy/d dz/d = 1. 1. 1. = relative spacing of dipoles in xtf, Ytf. ^tf directions 

• xox = TF coordinates XTp/d yrp/d ZTp/d corresponding to lattice site IXYZ= 

• (line containing comments) 

• dummy IXYZ {1 , 1) IXYZ(1,2) IXYZ(1,3) IC0MP(1,1) IC0MP(1,2) IC0MP(1,3) 

• dummy IXYZ {2, 1) IXYZ (2, 2) IXYZ (2, 3) IC0MP(2,1) ICOMP(2,2) ICOMP(2,3) 

• dummj/ IXYZ (3, 1) IXYZ(3,2) IXYZ(3,3) IC0MP(3,1) ICOMP(3,2) ICOMP(3,3) 



• dwmmj/ IXYZ (J, 1) IXYZ(J,2) IXYZ(J,3) IC0MP(J,1) IC0MP(J,2) IC0MP(J,3) 

• ... 

• dummy IXYZ (N, 1) IXYZ (N, 2) IXYZ (N, 3) ICOMP(N, 1) IC0MP(N,2) IC0MP(N,3) 

where dummy is a number (integer or floating point)that might, for example, identify the dipole. This 
number will not be used in any calculations. 

If the target material at location J is isotropic, ICOMP ( J, 1) , ICOMP (J, 2) , and ICOMP (J, 3) have 
the same value. 

demo file for target option FROM_FILE (homogeneous , isotropic target) 



8 








NAT 








1 , 


,000 







. 000 


0. 


,000 


= target vector al (in TF) 


0. 


,000 




1 


. 000 


. 


,000 


= target vector a2 (in TF) 


1 . 






1 




1 . 




= d_x/d d_y/d d_z/d (normally 


0. 


,5 




0, 


.5 


0. 


.5 


= X0(l-3) = location in lattice 


J 




JX 




JY 


JZ 


ICOMPX, ICOMPY, ICOMPZ 


1 















111 


2 












1 


1 : 


L 1 


3 









1 





1 : 


1 1 


4 









1 


1 


111 


5 




1 










111 


6 




1 







1 


111 


7 




1 




1 





111 


8 




1 




1 


1 


111 



The above sample target consists of 8 dipoles arranged to represent a cube. 

This example is homogeneous: All sites have composition 1 
The target origin XO is set to be at the center of the target 

Note that ICOMPX,ICOMPY,ICOMPZ could differ, allowing treatment of anisotropic targets, provided 
the dielectric tensor at each location is diagonal in the TF. 

demo file for target option FROM_FILE (inhomogeneous, isotropic target) — 

8 = NAT 

1.000 0.000 0.000 = target vector al (in TF) 

0. 000 1.000 0.000 = target vector a2 (in TF) 

1. 1. 1. = d_x/d d_y/d d_z/d (normally 111) 
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X0(l-3) = location in lattice of "target origin" 



This sample target consists of 8 dipoles arranged to represent a cube. 

This example is inhomogeneous: The lower half of the cube (JX=0) has isotropic composition 1 
The upper half of the cube (JX=1) has isotropic composition 2 
The target origin XO is set to be at the center of the target. 

Note that ICOMPX, ICOMPY, ICOMPZ can be different, allowing treatment of anisotropic targets, 
provided the dielectric tensor at each location is diagonal in the TF. 



20.1.1 Sample calculation in directory examples_exp/FROM_FILE 

Subdirectory examples_exp/FROM_FlLE contains ddscat . par for calculation of scattering by a 
0.5/im X Ifim x l//m Au block, represented by a 32 x 64 x 64 = 131072 dipole array, as well as the 
output files from the calculation. The target geometry is input via the file shape . dat. 

This target has V = 0.5/im^, and a^ff — (SV/At:)^/^ — 0.49237/im. The calculation is for an 
incident wavelength A = 0.50/im; the Au has refractive index m — 0.9656 + 1.8628i. The CCG method 
used is PBCGS2; the two orthogonal polarization require 29 and 30 iterations, respectively, to converge 
to the specified tolerance TOL = le-5. The computation used 165 MB of RAM, and required 208 
cpu sec on a 2.53 GHz cpu. 

N.B.: This is the same physical problem as the example in examples_exp/RCTGLPRSM (see 



§20.18.1 1, differing only in that in the present calculation the target geometry is input through the file 



shape . dat rather than generated by ddscat. 

20.2 ANIFRMFIL = General anistropic target defined by list of dipole loca- 
tions,"compositions", and material orientations obtained from a file 

This option causes DDSCAT to read the target geometry information from a file shape . dat instead 
of automatically generating one of the geometries listed below. The file shape . dat gives the num- 
ber N of dipoles in the target, the components of the "target axes" ai and k-i in the Target Frame 
(TF), the vector xo(l — 3) determining the correspondence between the integers IXYZ and actual co- 
ordinates in the TF, and specifications for the location and "composition" of each dipole. For each 
dipole J, the file shape . dat provides the location IXYZ (J, 1-3 ) , the composition identifier integer 
ICOMP (J, 1-3 ) specifying the dielectric function corresponding to the three principal axes of the di- 
electric tensor, and angles 9df, df, and /3df specifying the orientation of the local "Dielectric Frame" 



(DF) relative to the "Target Frame" (TF) (see 5 27 1. The DF is the reference frame in which the dielectric 
tensor is diagonalized. The Target Frame is the reference frame in which we specify the dipole locations. 

The shape.dat file is read by routine REASHP (file reashp.f90). The user can customize 
REASHP as needed to conform to the manner in which the target geometry is stored in file shape . dat. 
However, as supplied, REASHP expects the file shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements. 

• N = number of dipoles in target 

• Ola: o-iy o,iz = x,y,z Components (in Target Frame) of ai 

• 0'2x 0'2y o,2z = x,y,z Components (in Target Frame) of a2 
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• dx/d dy/d dz/d=\. 1. 1. = relative spacing of dipoles in xtf. Ytf. ztf directions 

• xqx xoy xoz = TF coordinates xtf / d utf /d ztf / d corresponding to lattice site IXYZ= 

• (line containing comments) 

• dummy IXYZ (1,1-3) ICOMP (1,1-3) THETADF(l) PHIDF(l) BETADF ( 1 ) 

• dummy IXYZ {2, 1-3) ICOMP (2, 1-3) THETADF(2) PHIDF(2) BETADF (2) 

• dummy IXYZ (3, 1-3) ICOMP (3, 1-3) THETADF(3) PHIDF(3) BETADF (3) 

• ... 

• dummj/ IXYZ (J, 1-3) ICOMP (J, 1-3) THETADF(J) PHIDF(J) BETADF (J) 

• ... 

• dwmmj/ IXYZ (N, 1-3) ICOMP (N, 1-3) THETADF (N) PHIDF (N) BETADF (N) 

Where dummy is a number (either integer or floating point) that might, for example, give the number 
identifying the dipole. This number will not be used in any calculations. 
THETADF PHIDF BETADF should be given in radians. 

Here is an example of the first few lines of a target description file suitable for target option AN I F RMF I L: 

demo file for target option ANIFRMFIL (this line is for comments) 
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This sample target consists of 8 dipoles arranged to represent a cube. 
Half of the cube (dipoles with JX=0) has isotropic composition 1 . For this case, the angles THETADF, 
PHIDF, BETADF do not matter, and it convenient to set them all to zero. 

The other half of the cube (dipoles with JX=1) consists of a uniaxial material, with dielectric function 2 
for E fields parallel to one axis (the "c-axis"), and dielectric function 3 for E fields perpendicular to the 
c-axis. The c-axis is 30" (0.5236 radians) away from xtf, and lies in the xtf-ztf plane (having been 
rotated by 1.5708 radians around xtf)- 

Note that 1C0MP(J,K) can be different for K=l,3, allowing treatment of anisotropic targets, provided 
the dielectric tensor at each location is diagonal in the TF. 



20.3 ANIELLIPS = Homogeneous, anisotropic ellipsoid 

SHPARi, SHPAR2, SHPAR3 define the ellipsoidal boundary: 



/ XTp/d 

[ SHPAR 



+ 



/ ZT^-p/d \ 1 
" 4 



VSHPAR2/ VSHPARs/ 
The TF origin is located at the centroid of the ellipsoid. 



(32) 
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20.4 ANI_ELL_2 = Two touching, homogeneous, anisotropic ellipsoids, with 
distinct compositions 

Geometry as for ELLIPS0_2; SHPARi, SHPAR2, SHPAR3 have same meanings as for ELLIPS0_2. 

Target axes ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf- 

Line connecting ellipsoid centers is || ai = xtf- 

TF origin is located between ellipsoids, at point of contact. 

It is assumed that (for both ellipsoids) the dielectric tensor is diagonal in the TF. User must set NC0MP=6 
and provide xx, yy, zz components of dielectric tensor for first ellipsoid, and xx, yy, zz components of 
dielectric tensor for second ellipsoid (ellipsoids are in order of increasing xtf)- 

20.5 ANI_ELL_3 = Three touching homogeneous, anisotropic ellipsoids with 
same size and orientation but distinct dielectric tensors 

SHPARi, SHPAR2, SHPAR3 have same meanings as for ELLIPS0_3. 

Target axis ai (1, 0, 0)tf (along fine of ellipsoid centers), and a2 = (0, 1, 0)tf- 

TF origin is located at center of middle elhpsoid. 

It is assumed that dielectric tensors are all diagonal in the TF. User must set NC0MP=9 and provide xx, 

yy, zz elements of dielectric tensor for first ellipsoid, xx, yy, zz elements for second elhpsoid, and xx, 
yy, zz elements for third elhpsoid (elhpsoids are in order of increasing jjtf)- 

20.6 ANIRCTNGL = Homogeneous, anisotropic, rectangular solid 

X, y, z lengths/o! = SHPARi, SHPAR2, SHPAR3. 

Target axes ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf in the TF. 

(a^TF, yxF, 2^tf) = (0, 0, 0) at middle of upper target surface, (where "up" = xtf)- (The target surface 

is taken to be d/2 about the upper dipole layer.) 

Dielectric tensor is assumed to be diagonal in the target frame. 

User must set NC0MP=3 and supply names of three files for e as a function of wavelength or energy: first 
for e^x, second for e^^, and third for e^^, as in the following sample ddscat . par file: 

' ============ Parameter file for v7 . 1 . ===================' 

PRELIMINARIES ****' 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGS2'= CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT'= CMETHD*6 (GPFAFT, FFTMKL) 
'GKDLDR'= CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 
' ANIRCTNGL' = CSHAPE*9 shape directive 

Initial Memory Allocation 
10 25 50 = upper bound on target extent 
' -k-k-k-k Target Geometry and Composition * + + + ' 
' ANIRCTNGL' = CSHAPE*9 shape directive 
10 25 50 = shape parameters SHPARI, SHPAR2 , SHPAR3 
3 = NCOMP = number of dielectric materials 

' /u/draine/DDA/diel/m2 . 00_0 . 10' = name of file containing dielectric function 
' /u/draine/DDA/diel/ml . 50_0 .00' 
' /u/draine/DDA/diel/ml . 50_0 .00' 

Error Tolerance ****' 
l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC| E>-ACA | X> ) / (NORM OF AC|E>) 

Interaction cutoff parameter for PEC calculations ****' 
5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 
' **** Angular resolution for calculation of <cos>, etc. ****' 
2.0 = ETASCA (number of angles is proportional to [ (2+x) /ETASCA] '^2 ) 

' **** Wavelengths (micron) ****' 

500. 500. 1 ' INV = wavelengths ( first, last, how many , how=LIN, INV, LOG) 

'k-k-k-k Effective Radii (micron) **** ' 

30.60 30.60 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 

'kkkk Define Incident Polarizations kkkk' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol. state) 

' Specify which output files to write ★★★★' 
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1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

Prescribe Target Rotations 
0. 0. 1 - BilAMI, BETAMX, NBETA (beta=rotation around al) 
0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
Specify first IWAV, IRAD, lORI (normally 0) ****' 
= first IWAV, first IRAD, first lORI (0 to begin fresh) 

Select Elements of S_ij Matrix to Print 
5 = NSMELTS = number of elements of S_ij to print (not more than 9) 

11 12 21 22 31 41 = indices ij of elements to print 

Specify Scattered Directions ****' 
'LFRAME' = CMDFRM*6 ('LFRAME' or 'TFRAME' for Lab Frame or Target Frame) 

2 = number of scattering planes 

0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 30 = phi, ... for plane B 

20.7 CONELLIPS = Two concentric ellipsoids 

SHPARi, SHPAR2, SHPAR3 = lengths/d of the outer ellipsoid along the xtf, Ytf, ztf axes; 
SHPAR4, SHPAR5, SHPARg = lengths/d of the inner ellipsoid along the xtf. Ytf, ztf axes. 
Target axes ai = (1, 0, 0)tf, a2 = (0, 1, 0)tf. 
TF origin is located at centroids of ellipsoids. 

The "core" within the inner ellipsoid is composed of isotropic material 1; the "mantle" between inner 
and outer ellipsoids is composed of isotropic material 2. 

User must set NC0MP=2 and provide dielectric functions for "core" and "mantie" materials. 

20.8 CYLINDERl = Homogeneous, isotropic finite cylinder 

SHPARi = length/d, SHPAR2 = diameter/d, with 

SHPAR3 = 1 for cyUnder axis ai || xtf: ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf; 
SHPAR3 = 2 for cylinder axis ai || ytf: ai = (0, 1, 0)tf and a2 = (0, 0, 1)tf; 
SHPAR3 = 3 for cylinder axis ai || ztf^ ^1 = (0, 0, 1)tf and a2 = (1, 0, 0)tf in the TF. 
TF origin is located at centroid of cylinder. 
User must set NC0MP=1. 

20.9 CYLNDRCAP = Homogeneous, isotropic finite cylinder with hemispheri- 
cal endcaps. 

SHPARi = cylinder length/d (not including end-caps!) and SHPAR2 = cylinder diameter/d, with cylinder 

axis = ai = (1,0,0)tf and a2 = (0, 1,0)tf- The total length along the target axis (including the 

endcaps) is (SHPARi+SHPAR2)d. 

TF origin is located at centroid of cylinder. 

User must set NC0MP=1. 

20.10 DSKRCTNGL = Disk on top of a homogeneous rectangular slab 

This option causes DDSCAT to create a target consisting of a disk of composition 1 resting on top of a 
rectangular block of composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic, 
ddscat . par should set NCOMP to 2 . 

The cylindrical disk has thickness SHPARi xd in the x-direction, and diameter SHPAR2X(i. The rect- 
angular block is assumed to have thickness SHPAR3xd in the x-direction, length SHPAR4X(i in the 
y-direction, and length SHPAR5 xd in the z-direction. The lower surface of the cylindrical disk is in the 
a; = plane. The upper surface of the slab is also in the a; = plane. 

The Target Frame origin (0,0,0) is located where the syrmnetry axis of the disk intersects the a; = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at {x/d,y/d, z/d) = {jx + 0.5, jy + Ay,jz + A^), where 
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jx, jy, and jz are integers. Aj, = or 0.5 depending on whether SHPAR4 is even or odd. = or 0.5 
depending on whether SHPAR5 is even or odd. 

Dipoles representing the disk are located at 
x/d = 0.5, 1.5, [int(SHPAR4 + 0.5) - 0.5] 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
Ooff = (3VT/47r)^/^, where Vr is the total volume of solid material in the target. For this geome- 
try, the number of dipoles in the target will be approximately N = [SHPARi x SHPAR2 x SHPAR3 + 
(7r/4) ((SHPAR4)^ x SHPAR5)], although the exact number may differ because of the dipoles are required 
to be located on a rectangular lattice. The dipole spacing d in physical units is determined from the 
specified value of acfi and the number N of dipoles in the target: d = {in/SNy^'^acff. This option 
requires 5 shape parameters: 
The pertinent line in dd scat. par should read 



SHPARI SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARi = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR2 = (disk diameter)/(i 

SHPAR3 = (brick thickness in xtf direction)/^ [material 2] 
SHPAR4 = (brick thickness in yxF direction)/^ 
SHPAR5 = (brick thickness in ztf direction)/o? 



The overall size of the target (in terms of numbers of dipoles) is determined by parameters ( SHPARI + SHPAR4 ) 
SHPAR2, and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
and SHPAR5. 

The physical size of the TUC is specified by the value of (in physical units, e.g. cm), specified in the 
file ddscat . par. 

The "computational volume" is determined by (SHPARi + SHPAR4)x SHPAR2 x SHPAR3. 
The target axes (in the TF) are set to ai = xtf = (17 0,0)tf - i-e., normal to the "slab" - and 
^2 — Ytf — (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all other 



targets, set by the usual orientation angles Q, $, and (3 (see 518 above); for example, 6 = would be 
for radiation incident normal to the slab. 



20.11 DW1996TAR = 13 block target used by [Draine & Weingartner] ( |1996l ). 

Single, isotropic material. Target geometry was used in study by |Draine & Weingartner| ( fl996| l of radia- 
tive torques on irregular grains, ai and a2 are principal axes with largest and second-largest moments of 
inertia. User must set NC0MP = 1. Target size is controlled by shape parameter SHPAR ( 1 ) = width of 
one block in lattice units. 
TF origin is located at centroid of target. 



20.12 ELLIPSOID = Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARi, SHPAR2, SHPAR3 in the x, y, z directions in the TF: 

( ^Z^V + ( V+ ( ^I^V = i (33 

VSHPARirfy VSHPAR2d/ VSHPAR3d/ 4 ' 

where d is the interdipole spacing. 

The target axes are set to ai = (1, 0, 0)tf and 3.2 — (0, 1, 0)tf- 

Target Frame origin = centroid of ellipsoid. 

User must set NC0MP=1 on line 9 of ddscat . par. 

A homogeneous, isotropic sphere is obtained by setting SHPARi = SHPAR2 = SHPAR3 = diameter/d. 
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Figure 8: Normalized electric field intensity |E|^/|Eo|^ along a line parallel to the direction of propagation, and passing 
through the center of an Au sphere of radius a = 0.3979/im, for light with wavelength A = 0.5/im. The calculation in exam- 
ples_exp/ELLIPSOID_NEARFIELD was done with dipole spacing d = a/48.49 = 0.00821/im. 



20.12.1 Sample calculation in directory examples_exp/ELLIPSOID 

The directory examples_exp/ELLIPSOID contains ddscat .par for calculation of scattering by 
a sphere with refractive index m — 1.5 + O.Oli and 27ra/A = 5, represented by a = 59728 dipole 
pseudosphere just fitting within a 48 x 48 x 48 computational volume, as well as the output files from 
the calculation. The calculation with 2na/X — 5 has \m\kd = 0.309. The computation used 144 MB of 
RAM and required 63 cpu sec on a 2.53 GHz cpu. 

20.12.2 Sample calculation in directory examples_exp/ELLIPSOID_NEARFIELD 

The directory examples_exp/ELLIPSOID_NEARFIELD contains ddscat . par for calculation of 
(1) far-field scattering and (2) E in and near the target for a sphere with refractive index m — 0.96+l.Oli 
(refractive index of Au at A = 0.5/im) and 27ra/A = 5. The spherical target is represented by a 
N = 59728 dipole pseudosphere just fitting within a 48 x 48 x 48 computational volume, as well as the 
output files from the calculation. 

In physical units with A = 0.5/im, Ooff = 5 x A/27r = 0.39789/j,m. The calculation with 27ra/A = 5 
has |m|fc(i = 0.309. The computation used 144 MB of RAM and required 60 cpu sec on a 2.53 GHz 
cpu. 

The nearfield calculation is specified to extend throughout a computational volume extending the 
original A8d x 48d x ASd computational volume by 50% in all directions, to become a 96d x 96d x 96d 
volume centered on the sphere. E is evaluated at all points in this volume. The nearfield calculation 
used 62 MB of RAM and required just 9.6 cpu sec. The nearfield calculation creates the binary files 
w000r000k000.Elandw000r000k000.E2,one for each of the two incident polarizations. 

After the nearfield calculation is complete, the program readnf is used to read the file wOOOrOOOkOOO.El 
(specified in readnf . par) and extract E at 501 points along a line specified in readnf . par - the 
line runs along the xtf axis through the center of the sphere. The results are shown in Figure |8] 

20.13 ELLIPSO_2 = Two touching, homogeneous, isotropic eUipsoids, with dis- 
tinct compositions 



SHPARi, SHPAR2, SHPAR3=x-length/rf, y-length/d, z-length/d of one ellipsoid. The two ellipsoids 
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have identical shape, size, and orientation, but distinct dielectric functions. The line connecting ellipsoid 
centers is along the xxF-axis. Target axes ai = (1,0, 0)tf [along line connecting ellipsoids] and 

a2 = (0,l,0)TF. 

Target Frame origin = midpoint between ellipsoids (where ellipsoids touch). 

User must set NC0MP=2 and provide dielectric function file names for both elhpsoids. Ellipsoids are 
in order of increasing xtf- first dielectric function is for elUpsoid with center at negative xtf, second 
dielectric function for ellipsoid with center at positive a;TF- 

20.14 ELLIPSO_3 = Three touching homogeneous, isotropic ellipsoids of equal 
size and orientation, but distinct compositions 

SHPARi, SHPAR2, SHPAR3 have same meaning as for ELLIPS0_2. Line connecting ellipsoid cen- 
ters is parallel to xtf axis. Target axis ai = (1,0,0)tf (along Une of ellipsoid centers), and a2 = 

(0,1,0)tf. 

Target Frame origin = centroid of middle ellipsoid. 

User must set NC0MP=3 and provide (isotropic) dielectric functions for first, second, and third ellipsoid. 

20.15 HEX_PRISM = Homogeneous, isotropic hexagonal prism 

SHPARi = (Length of prism)/d = (distance between hexagonal faces)/d, 

SHPAR2 = (distance between opposite vertices of one hexagonal face)/d = 2xhexagon side/d. 

SHPAR3 selects one of 6 orientations of the prism in the Target Frame (TF). 

Target axis ai is along the prism axis (i.e., normal to the hexagonal faces), and target axis a2 is normal 
to one of the rectangular faces. There are 6 options for SHPAR3: 

SHPAR3 = 1 for ai II Xtf and a2 || yxF ; SHPAR3 = 2 for ai || xtf and a2 || ztf ; 
SHPAR3 = 3 for ai II ytf and 3.2 || xtf ; SHPAR3 = 4 for ai || yxF and a2 || ztf ; 
SHPAR3 = 5 for ai II Ztf and a2 || xtf ; SHPAR3 = 6 for ai || ztf and a2 || yxF 
TF origin is located at the centroid of the target. 
User must set NC0MP=1. 

20.16 LAYRDSLAB = Multilayer rectangular slab 

Multilayer rectangular slab with overall x, y, z lengths = SHPARi x d 

tty = SHPAR2 X d, 

= SHPAR3 X d. 

Upper surface is at xtf = 0, lower surface at xtf = —SHPARi x d 

SHPAR4 = fraction which is composition 1 (top layer). 

SHPAR5 = fraction which is composition 2 (layer below top) 

SHPARe = fraction which is composition 3 (layer below comp 2) 

1 - (SHPAR4 + SHPAR5 + SHPARg) = fraction which is composition 4 (bottom layer). 

To create a bilayer slab, just set SHPAR5 = SHPARe = 

To create a trilayer slab, just set SHPARg = 

User must set NCOMP=2,3, or 4 and provide dielectric function files for each of the two layers. Top 
dipole layer is at xtf = —d/2. Origin of TF is at center of top surface. 

20.17 MLTBLOCKS = Homogeneous target constructed from cubic "blocks" 

Number and location of blocks are specified in separate file blocks . par with following structure: 
one line of comments (may be blank) 
PRIN (= or 1 - see below) 
N (= number of blocks) 

B (= width/rf of one block) 

a;TF TJTF Ztf (= position of 1st block in units of Bd) 
Xtf Vtf Ztf (= position of 2nd block in units of BcQ ) 



TARGET GENERATION: ISOLATED FINITE TARGETS 



Xtf Utf -2tf (= position of Nth block in units of Bd) 
If PRIN=0, then ai = (1, 0, 0)tf, a2 = (0, 1, 0)tf- If PRIN=1, then ai and a2 are set to principal axes 
with largest and second largest moments of inertia, assuming target to be of uniform density. User must 
setNCOMP=l. 



20.18 RCTGLPRSM = Homogeneous, isotropic, rectangular solid 

X, y, z lengths/d = SHPARi, SHPAR2, SHPAR3. 
Target axes ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf- 

TF origin at center of upper surface of solid: target extends from xtf/cI = — SHPARi to 0, 
VTF/d from -0.5 x SHPAR2 to +0.5 x SHPARj 
^;TF/rffrom -0.5 x SHPAR3 to +0.5 x SHPAR3 
User must set NC0MP=1. 



20.18.1 Sample calculation in directory examples_exp/RCTGLPRSM 

The directory examples_exp/RCTGLPRSM contains ddscat . par for calculation of scattering by 
a 0.5/im x 1/im x l/im Au block, represented by a 32 x 64 x 64 dipole array, together with output files 
from the calculation. The Au has refractive index m = 0.9656 + 1.8628*. The DDA calculation has 
\m\kd = 0.4120. The calculation used 165 MB of RAM, and required 207 cpu sec on a 2.53 GHz cpu. 



20.19 
axis. 



RCTGLBLK3 = Stack of 3 rectangular blocks, with centers on the xtf 



Each block consists of a distinct material. There are 9 shape parameters: 

SHPARi = (upper solid thickness in xtf direction)/^ [material 1] 

SHPAR2 = (upper solid width in yxF direction)/^ 

SHPAR3 = (upper solid width in Ztf direction)/(i 

SHPAR4 = (middle solid thickness in xtf direction)/d [material 2] 

SHPAR5 = (middle solid width in yxp direction)/^ 

SHPARe = (middle sohd width in ztf direction)/(i 

SHPAR7 = (lower solid thickness in xtf direction)/^ [material 3] 

SHPARg = (lower solid width in yxF direction)/^ 

SHPARg = (lower solid width in ztf direction)/(i 

TF origin is at center of top surface of material 1. 



20.20 SLAB_HOLE = Rectangular slab with a cylindrical hole. 

The target consists of a rectangular block with a cylindrical hole with the axis passing through the 
centroid and aligned with the xtf axis. The block dimensions are a x 6 x c. The cyhndrical hole has 
radius r. The pertinent line in ddscat . par should read 
SHPARi SHPAR2 SHPAR3 SHPAR4 
where 

SHPARi = a/did is the interdipole spacing) 

SHPAR2 = b/a 
SHPAR3 = c/a 
SHPAR4 = r/a 

Ideally, SHPARi, SEP AR2 x SHPARi, SHPAR3 x SHPARi will be integers (so that the cubic lattice can 
accurately approximate the desired target geometry), and SHPAR4X SHPARi will be large enough for 
the circular cross section to be well-approximated. 

The TF origin is at the center of the top surface (the top surface lies in the yxF ^ ^tf plane, and 
extends from j/tf = —b/2 to +6/2, and ztf = — c/2 to +c/2). The cylindrical hole axis runs from 
(a;TF = 0, Utf = 0, ztf = 0) to {xtf = -a, utf = 0, ztf = 0). 
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20.21 SPHERES_N = Multisphere target = union of N spheres of single isotropic 
material 

Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an exter- 
nal file, whose name (enclosed in single quotes) is passed through ddscat . par. The length of the file 
name should not exceed 80 characters. The pertinent line in ddscat . par should read 

SHPARi SHPAR2 'filename' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR2= to have oi = (1, 0, 0)tf, ^2 = (0, 1, 0)tf. 

SHPAR2= 1 to use principal axes of moment of inertia tensor for ai (largest /) and a2 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPARi, which is the extent of the multisphere target in the a;-direction, in units of the lattice spacing d. 
The file filename' should have the following structure: 

N (= number of spheres) 
line of comments (may be blank) 
line of comments (may be blank) [N.B 
line of comments (may be blank) [N.B 
line of comments (may be blank) [N.B 
xi yi z\ a\ (arb. units) 
xi 2/2 -22 a2 (arb. units) 

yw zn (In (arb. units) 

where Xj, yj, zj are the coordinates (in the TF) of the center of sphere j, and aj is the radius of sphere j. 
Note that Xj, yj, Zj, aj (J = 1, N) establish only the shape of the A?^— sphere target. For instance, 
a target consisting of two touching spheres with the line between centers parallel to the x axis could 
equally well be described by lines 6 and 7 being 

0.5 

1 0.5 

or 

1 

2 1 

The actual size (in physical units) is set by the value of agg specified in ddscat . par, where, as always, 
aeff = (SV/Att)^/^, where V is the total volume of material in the target. 

Target axes ai and 3.2 are set to be principal axes of moment of inertia tensor (for uniform density), 
where ai corresponds to the largest eigenvalue, and a2 to the intermediate eigenvalue. 
The TF origin is taken to be located at the volume-weighted centroid. 
User must set NC0MP=1. 

20.22 SPHROID_2 = Two touching homogeneous, isotropic spheroids, with dis- 
tinct compositions 

First spheroid has length SHPARi along symmetry axis, diameter SHPAR2 perpendicular to symmetry 
axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 perpendicular to 
symmetry axis. Contact point is on line connecting centroids. Line connecting centroids is in Xtf 
direction. Symmetry axis of first spheroid is in ytf direction. Synmietry axis of second spheroid is in 
direction yxp cos(SHPAR5) + ztf sin(SHPAR5), and SHPAR5 is in degrees. If SHPARg = 0., then target 
axes ai = (1, 0, 0)tf, ^2 = (0, 1, 0)tf- If SHPARg = 1., then axes ai and a2 are set to principal axes 
with largest and 2nd largest moments of inertia assuming spheroids to be of uniform density. 
Origin of TF is located between spheroids, at point of contact. 
User must set NC0MP=2 and provide dielectric function files for each spheroid. 



.: changed from v7. 0.7] 
.: changed from v7.0.7] 
.: changed from v7.0.7] 
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20.23 SPH_ANI_N = Multisphere target consisting of the union of spheres of 
various materials, possibly anisotropic 

Spheres may NOT overlap. The relative locations and sizes of these spheres are defined in an external 
file, whose name (enclosed in single quotes) is passed through ddscat .par. The length of the file 
name should not exceed 80 characters. Target axes ai and a.2 are set to be principal axes of moment 
of inertia tensor (for uniform density), where ai corresponds to the largest eigenvalue, and 3.2 to the 
intermediate eigenvalue. 

The TF origin is taken to be located at the volume-weighted centroid. 
The pertinent line in ddscat . par should read 
SHPARi SHPAR2 'filename' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR2= to have oi = (1, 0, 0)tf, 02 = (0, 1, 0)tf in Target Frame. 

SHPAR2= 1 to use principal axes of moment of inertia tensor for ai (largest /) and a2 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPARi, which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. 
The file filename' should have the following structure: 



TV (= number of spheres) 

line of comments (may be blank) 

line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 

xi ui zi ri Cxi Cyi Czi 0df,i ^'df,i /3df,i 
X2 V2 Z2 r2 Cx2 Cy2 Cz2 6df,2 *df,2 /3df,2 

xn UN zn tn Cxm CyN Czn Qdf.n ^df,n /3df,n 



where Xj, yj, Zj are the coordinates of the center, and rj is the radius of sphere j (arbitrary units), 
Cxj, Cyj, Czj are integers specifying the "composition" of sphere j in the x,y,z directions in the 
"Dielectric Frame" (see 5 27 1 of sphere j, and &uF.j ^df j /^dfj are angles (in radians) specifying 
orientation of the dielectric frame (DF) of sphere j relative to the Target Frame. Note that Xj, yj, zj, rj 
{j = 1, N) establish only the shape of the A^— sphere target, just as for target option NSPHER. The 
actual size (in physical units) is set by the value of a^ff specified in ddscat .par, where, as always, 
Ocff = (3y/47r)^/'^, where V is the volume of material in the target. 

User must set NCOMP to the number of different dielectric functions being invoked (i.e., the range of 

{Cxj , Cyj , Czj }. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file 'filename ' would be 
simply 



TV (= number of spheres) 

line of comments (may be blank) 

line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 
line of comments (may be blank) [N.B.: changed from v7.0.7] 

2/1 zi ri 1 1 1 
a;2 2/2 ^2 ^2 1 1 1 



Xn i/AT z^v rjv 1 1 1 
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20.23.1 Sample calculation in directory examples_exp/SPH_ANI_N 

Subdirectory examples_exp/SPH_ANI_N contains ddscat.par for calculating scattering by a 
random aggregate of 64 spheres [aggregated according following the "BAM2" aggregation process de- 
scribed by ^Shen et al., ( ,2008 ) ]. 32 of the spheres are assumed to consist of "astrosilicate", and 32 of 
crystalline graphite, with random orientations for each of the 32 graphite spheres. Each sphere is as- 
sumed to have a radius 0.050/im. The entire cluster is represented by iV = 7947 dipoles, or about 124 
dipoles per sphere. Scattering and absorption are calculated for A — 0.55/im. 

20.24 TETRAHDRN = Homogeneous, isotropic tetrahedron 

SHPARi=length/d of one edge. Orientation: one face parallel to yxF, ztf plane, opposite "vertex" is 
in +XTF direction, and one edge is parallel to ztf. Target axes ai = (1, 0, 0)tf [emerging from one 
vertex] and 3.2 = (0, 1, 0)tf [emerging from an edge] in the TF. User must set NC0MP=1. 

20.25 TRNGLPRSM = Triangular prism of homogeneous, isotropic material 

SHPARi, SHPAR2, SHPAR3, SEP AR4 = a/d, b/a, c/a, L/a 

The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the lattice 
spacing. The triangular cross-section has interior angles a, (3, 7 (opposite sides a, b, c) given by cos a — 
(fe^ + — a?)/2bc, cos/3 = (a^ + — 6^)/2ac, cos 7 = (a^ + fe^ — c^)/2ab. In the Target Frame, 
the prism axis is in the x direction, the normal to the rectangular face of width a is (0,1,0), the normal to 
the rectangular face of width b is (0, — cos 7, sin 7), and the normal to the rectangular face of width c is 
(0, — cos/3, — sin/3). 

20.26 UNIAXICYL = Homogeneous finite cylinder with uniaxial anisotropic di- 
electric tensor 

SHPARi, SHPAR2 have same meaning as for CYLINDERl. Cylinder axis = ai = (1,0, 0)tf, a2 = 
(0, 1, 0)tf- It is assumed that the dielectric tensor e is diagonal in the TF, with eyy — e^z- User must set 
NC0MP=2. Dielectric function 1 is for E || ai (cylinder axis), dielectric function 2 is for E _L ai. 



20.27 Modifying Existing Routines or Writing New Ones 

The user should be able to easily modify these routines, or write new routines, to generate targets with 
other geometries. The user should first examine the routine target . f 90 and modify it to call any 
new target generation routines desired. Alternatively, targets may be generated separately, and the target 
description (locations of dipoles and "composition" corresponding to x,y,z dielectric properties at each 
dipole site) read in from a file by invoking the option FROM_FILE in ddscat . f 90. 

Note that it will also be necessary to modify the routine reapar.f90so that it will accept whatever 
new target option is added to the target generation code . 



20.28 Testing Target Generation using CALLTARGET 

It is often desirable to be able to run the target generation routines without running the entire DDSCAT 
code. We have therefore provided a program CALLTARGET which allows the user to generate targets 
interactively; to create this executable just typ^***] 

Non-Linux sites: The source code for CALLTARGET is in the file CALLTARGET . f 90. You must compile and link 

CALLTARGET . f 90, ddcommon . f 90, dsYevj3.f90, errmsg.f90, gasdev.f90, p_lm.f90, prinaxis . f 90, 
ran3.f90, reashp.f90, sizer.f90, tar2el.f90, tar2sp.f90, tar3el.f90, taranirec . f 90, 
tarblocks . f 90, tarcel.f90, tarcyl.fgo, tarcylcap . f 90, tarell.f90, target. f90, targspher . f 90, 
tarhex.f90, tarnas.f90, tarnsp.f90, tarpbxn.f90, tarprsm. f 90, tarrctblkS . f 90, tarrecrec . f 90, 
tarslblin. f 90, tartet . f 90, andwrimsg. f 90. 
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make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may need 
to edit the code to change the device number IDVOUT as for DDSCAT (see ^6.7| above). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target . out is the same as the format 
of the shape . dat files read if option FROM_FILE is used (see above). Therefore you can simply 

mv target. out shape.dat 
and then use DDSCAT with the option FROM_FILE (or option ANIFRMFIL in the case of anisotropic 
target materials with arbitrary orientation relative to the Target Frame) in order to input a target shape 
generated by CALLTARGET. 

Note that CALLTARGET - designed to generate finite targets - can be used with some of the "PBC" 
target options (see ^2T| below) to generate a list of dipoles in the TUC. At the moment, CALLTARGET 
has support for target options BISLINPBC, DSKBLYPBC, and DSKRCTPBC. 



21 Target Generation: Periodic Targets 



A periodic target consists of a "Target Unit Cell" (TUC) which is then repeated in either the ytf direc- 
tion, the ztf direction, or both. Please see Draine & Flatau ( 2008| l for illustration of how periodic targets 
are assembled out of TUCs, and how the scattering from these targets in different diffraction orders M 
or (M, A^) is constrained by the periodicity. 

The following options for the TUC geometry are included in DDSCAT: 



FRMFILPBC : TUC geometry read from file (5 21 . 1 



ANIFILPBC : TUC geometry read from file, anisotropic materials supported ( ^21.2 1 



BISLINPBC : TUC = bilayer slab (§21.3 1 



CYLNDRPBC : TUC = finite cylinder (^2L4i 
DSKBLYPBC : TUC = disk plus bilayer slab ( pT3 1 



DSKRCTPBC : TUC = disk plus brick (^ITS i 
HEXGONPBC : TUC = hexagonal prism (pLT]) 



LYRSLBPBC : TUC = layered slab (up to 4 layers) ( ^218 1 
RCTGL_PBC : TUC = brick 

RECRECPBC : TUC = brick resting on brick ( ^200] ) 
SLBHOLPBC : TUC = brick with cylindrical hole ( pTTT) 



SPHRN_PBC : TUC = N spheres ( pLlI i 
TRILYRPBC : TUC = three stacked bricks (^JTOjl 
Each option is decribed in detail below. 



21.1 FRMFILPBC = periodic target with TUC geometry and composition input 
from a file 

The TUC can have arbitrary geometry and inhomogeneous composition, and is assumed to repeat peri- 
odically in either 1-d (y or z) or 2-d (y and z). 
The pertinent line in ddscat . par should read 
SHPARi SHPAR2 'filename' (quotes must be used) 
SHPARi = Ptj/d (Py = periodicity in yxF direction) 
SHPAR2 = Pz/d (Pz = periodicity in ztf direction) 

filename is the name of the file specifying the locations of the dipoles, and the "composition" at each 
dipole location. The composition can be anisotropic, but the dielectric tensor must be diagonal in the 
TP. The shape and composition of the TUC are provided exactly as for target option FROM_FILE - see 
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If SHPARi = then the target does not repeat in the yxF direction. 
If SHPAR2 = then the target does not repeat in the ztf direction. 

21.2 ANIFILPBC = general anisotropic periodic target with TUC geometry and 
composition input from a file 

The TUC can have arbitrary geometry and inhomogeneous composition, and is assumed to repeat peri- 
odically in either 1-d (y or z) or 2-d (y and z). 
The pertinent line in dd scat. par should read 
SHPARi SHPAR2 'filename' (quotes must be used) 
SHPARi = Py/d (Py = periodicity in yxF direction) 
SHPAR2 = Pz/d (Pz = periodicity in ztf direction) 

filename is the name of the file specifying the locations of the dipoles, and the "composition" at each 
dipole location. The composition can be anisotropic, and the dielectric tensor need not be diagonal in 
the TF. The shape and composition of the TUC are provided exactly as for target option ANIFRMFIL - 
see ^02] 

If SHPARi = then the target does not repeat in the yxF direction. 
If SHPAR2 = then the target does not repeat in the ztf direction. 



21.3 BISLINPBC = Bi-Layer Slab with Parallel Lines 

The target consists of a bi-layer slab, on top of which there is a "line" with rectangular cross-section. 
The "line" on top is composed of material 1, has height Xi (in the Xtf direction), width Yi (in the yxF 
direction), and is infinite in extent in the ztf direction. 

The bilayer slab has width Y2 (in the yxF direction). It is consists of a layer of thickness X2 of material 
2, on top of a layer of material 3 with thickness X^. 



SHPARi 


= Xr/d 


(Xi 


= thickness of line) 


SHPAR2 


= Y,/d 


(Yi = 


: width of line) 


SHPAR3 


= X2/d 


(X2 


= thickness of upper layer of slab) 


SHPAR4 


= X^/d 


(Xs 


= thickness of lower layer of slab) 


SHPAR5 


= Y2/d 


(Y2 = 


: width of slab) 


SHPARe 


= Py/d 


(Py-- 


= periodicity in yxF direction). 



If SHPARg = 0, the target is NOT periodic in the yxF direction, consisting of a single column, infinite 
in the ztf direction. 

21.4 CYLNDRPBC = Target consisting of homogeneous cylinder repeated in 
target y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an infinite array of cylinders. The individual 
cylinders are assumed to be homogeneous and isotropic, just as for option RCTNGL (see ^20.18| l. 
Let us refer to a single cylinder as the Target Unit Cell (TUC). The TUC is then repeated in the target y- 
and/or z-directions, with periodicities PYDxd and PZDxd, where d is the lattice spacing. To repeat in 
only one direction, set either PYD or PZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 

where SHPARi, SHPAR2, SHPAR3, SHPAR4, SHPAR5 are numbers: 
SHPARi = cylinder length along axis (in units of d) in units of d 
SHPAR2 = cylinder diameter/c? 
SHPAR3 = 1 for cylinder axis || xtf 
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= 2 for cylinder axis || yxF 

= 3 for cylinder axis || ztf (see below) 
SHPAR4 = PYD = periodicity/d in yxF direction ( = to suppress repetition) 
SHPAR5 = PZD = periodicity/d in ztf direction ( = to suppress repetition) 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARl 
and SHPAR2. The orientation of a single cylinder is determined by SHPAR3. The periodicity in the TF 
y and z directions is determined by parameters SHPAR4 and SHPAR5. 

The physical size of the TUC is specified by the value of a^s (in physical units, e.g. cm), specified in 
the file ddscat .par, with the usual correspondence d — {in/SNy/^acs, where N is the number of 
dipoles in the TUC. 

With target option CYLNDRPBC, the target becomes a periodic structure, of infinite extent. 

• If NPY > and NPZ > 0, then the target cylindrical TUC repeats in the yxF direction, with 
periodicity NPY x d. 

• If NPY = and NPZ > then the target cylindrical TUC repeats in the Ztf direction, with 
periodicity NPZ x d. 

• If NPY > and NPZ > then the target cylindrical TUC repeats in the yxF direction, with 
periodicity NPY x d, and in the ztf direction, with periodicity NPZ x d. 

Target Orientation: The target axes (in the TF) are set to 3.1 — (1, 0, 0)tf and a2 — (Oi I7 0)tf- Note 
that ai does not necessarily coincide with the cylinder axis: individual cylinders may have any of 3 
different orientations in the TF. 



Example 1: One could construct a single infinite cylinder with the following two lines in dds cat .par: 



100 1 100 

1.0 100.49 2 1.0 0. 



The first line ensures that there will be enough memory allocated to generate the target. The TUC 
would be a thin circular "slice" containing just one layer of dipoles. The diameter of the circular sUce 
would be about 100.49d in extent, so the TUC would have approximately (7r/4) x (100.49)^ = 7931 
dipoles (7932 in the actual realization) within a 100 x 1 x 100 "extended target volume". The TUC would 
be oriented with the cylinder axis in the yxF direction (SHPAR3=2) and the structure would repeat in 
the yxF direction with a period of 1.0 x d. SHPAR5=0 means that there will be no repetition in the z 
direction. As noted above, the "target axis" vector ai = xtf- 

Note that SHPARi, SHPAR2, SHPAR4, and SHPAR5 need not be integers. However, SHPAR3, deter- 
mining the orientation of the cylinders in the TF, can only take on the values 1,2,3. 

The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles /3, 0, and $ (see § 18 above); for example, 9 = would be for radiation incident 
normal to the periodic structure. 



21.4.1 Sample calculation in directory examples_exp/CYLNDRPBC 

The subdirectory examples_exp/CYLNDRPBC contains ddscat . par for calculating scattering by 
an infinite cylinder with m = 1.33 + O.Olz for 27tR/X = 10, where R is the cylinder radius. The TUC is 
a disk of thickness d. The sample calculation calls for the cylinder diameter 2R to be 64.499(i, where d 
is the dipole spacing, so that the computational volume With this choice, the TUC (disk of thickness d) 
turns out to be represented by iV = 3260 dipoles, and the DDA calculation is done with \m\kd ~ 0.4129. 

With the standard error tolerance T0L=1 . Oe-5, the calculation converges in IT=16 iterations for 
one incident polarization, I T= 1 3 iterations for the orthogonal polarization. The scattering properties are 
reported inwOOOrOOOkOOO.sca. The calculation used required 40 MB of RAM, and used 39 cpu sec 
on a 2.53 GHz cpu. 
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21.5 DSKBLYPBC = Target consisting of a periodic array of disks on top of a 
two-layer rectangular slabs. 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of Target Unit 
Cells (TUCs), each TUC consisting of a disk of composition 1 resting on top of a rectangular block con- 
sisting of two layers: composition 2 on top and composition 3 below. Materials 1, 2, and 3 are assumed 
to be homogeneous and isotropic. 
This option requires 8 shape parameters: 
The pertinent line in dds cat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARg SHPAR7 SHPARg 

where 

SHPARi = disk thickness in x direction (in Target Frame) in units of d 

SHPAR2 = (disk diameter)/d 

SHPAR3 = (upper slab thickness)/(i 

SHPAR4 = (lower slab thickness)/^ 

SHPAR5 = (slab extent in yxF direction)/d 

SHPARg = (slab extent in ztf direction)/rf 

SHPAR7 = period in yxF direction/d 

SHPARg = period in ztf direction/d 

The physical size of the TUC is specified by the value of (in physical units, e.g. cm), specified 
in the file dds cat . par. 

The "computational volume" is determined by 
(SHPARi + SHPAR3 + SHPAR4)x SHPAR5 X SHPARg. 

The lower surface of the cylindrical disk is in the a; = plane. The upper surface of the slab is also 
in the a; = plane. It is required that SHPAR2 < min(SHPAR4, SHPAR5). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the a; = 
plane (the upper surface of the slab, and the lower surface of the disk). 
In the Target Frame, dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPARi + 0.5) - 0.5] 
and at (y, z) values 

y/d = ±0.5,±1.5,...and 

z/d = ±0.5, ±1.5, ... satisfying 

(y2 + ^2) < (sHPAR2/2)2d2. 

Dipoles representing the rectangular slab are located at (x/d, y/d, z/d) — [jx + 0.5, jy + Ay, jz + 
Az), where jx, jy, and jz are integers. Ay = or 0.5 depending on whether SHPAR5 is even or odd. 
A.Z = or 0.5 depending on whether SHPARg is even or odd. 

The TUC is repeated in the target y- and z-directions, with periodicities SHPARyxd and SHPARgxd. 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
aeff = (3Vtuc / 47r) ^Z'', where Vrijc is the total volume of solid material in one TUC. For this geometry, 
the number of dipoles in the target will be approximately 

N = (77/4) X SHPARi X (SHPARa)^ + [SHPAR3 + SHPAR4] x SHPAR5 x SHPARg 

although the exact number may differ because of the dipoles are required to be located on a rectangular 
lattice. The dipole spacing d in physical units is determined from the specified value of a^s and the 
number N of dipoles in the target: d = (47r/3A'')^/^aefi. 

The target axes (in the TF) are set to ai = Xtf = (Ij 0, 0)tf - i-c, normal to the "slab" - and 3.2 = 
yxF = (0, 1, 0)tf. The orientation of the incident radiation relative to the target is, as for all other 
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targets, set by the usual orientation angles /3, O, and $ (see 518 above); for example, 8 = would be 
for radiation incident normal to the slab. 

21.6 DSKRCTPBC = Target consisting of homogeneous rectangular brick plus 
a disk, extended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of Target Unit Cells. 
Each Target Unit Cell (TUC) consists of a disk of composition 1 resting on top of a rectangular block of 
composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic. 

The cylindrical disk has thickness SHPARixd in the x-direction, and diameter SHPAR2X(i. The rect- 
angular block is assumed to have thickness SHPARsxd in the x-direction, extent SHPAR4X(i in the 
y-direction, and extent SHPARsxd in the z-direction. The lower surface of the cylindrical disk is 
in the x = plane. The upper surface of the slab is also in the a; = plane. It is required that 
SHPAR2<min(SHPAR4, SHPAR5). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the x = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at {x/d, y/d, z/d) — {jx + 0.5, jj, + Ay,jz + Az), where 
jx, jy, and jz are integers. Ay = or 0.5 depending on whether SHPAR4 is even or odd. Az — or 0.5 
depending on whether SHPAR5 is even or odd. 
jx = -[int(SHPAR3 + 0.5)],...,-l. 

jy = -[int(0.5 X SHPAR4 - 0.5) + 1] , ... , = -[int(SHPAR4 + 0.5) - 0.5] 

y/d = -[int(0.5 x SHPAR4 + 0.5) - 0.5], [int(0.5 x SHPAR4 + 0.5) - 0.5] 

z/d = -[int(0.5 x SHPAR5 + 0.5) - 0.5], [int(0.5 x SHPAR5 + 0.5) - 0.5] 
where int(a;) is the greatest integer less than or equal to x. Dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPAR4 + 0.5) - 0.5] 
and at {y, z) values 

2//d= ±0.5,±1.5,... and 

z/d — ±0.5, ±1.5, ... satisfying 

(2/^ + z^) < (SHPAR5/2)2(i2. 
The TUC is repeated in the target y- and z-directions, with periodicities SHPARgxd and SHPARyxd. 
As always, the physical size of the target is fixed by specifying the value of the effective radius Ooff = 
(3VTUc/47r)^/'^, where Vtuc is the total volume of solid material in one TUC. For this geometry, 
the number of dipoles in the target will be approximately N = [SHPARi x SHPAR2 x SHPAR3 + 
(7r/4)((SHPAR4)^ x SHPAR5)], although the exact number may differ because of the dipoles are re- 
quired to be located on a rectangular lattice. The dipole spacing d in physical units is determined from 
the specified value of aoff and the number N of dipoles in the target: d — (An/SNy^'^acs. This option 
requires 7 shape parameters: 
The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARe SHPAR7 
where 

SHPARi = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR2 = (disk diameter)/^ 

SHPAR3 = (brick thickness in xtf direction)/^ [material 2] 
SHPAR4 = (brick length in yTF direction)/^ 
SHPAR5 = (brick length in ztf direction)/^ 
SHPARe = periodicity in yxF direction/d 
SHPAR7 = periodicity in zxf direction/d 

The overall extent of the TUC (the "computational volume") is determined by parameters ( SHPARI +SHPAR4 ) 
max ( SHPAR2 , SHPAR4, and max ( SHPAR3 , SHPAR5 ) . The periodicity in the TF y and z directions is 
determined by parameters SHPARg and SHPAR7. 

The physical size of the TUC is specified by the value of OeS (in physical units, e.g. cm), specified in the 
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file ddscat . par. 

The target is a periodic structure, of infinite extent in the target y- and z- directions. The target axes 
(in the TF) are set to ai = xxf = (Ij 0, 0)tf - i-e-. normal to the "slab" - and a2 — yxF = (0, 1, 0)tf- 
The orientation of the incident radiation relative to the target is, as for all other targets, set by the usual 
orientation angles /3, 8, and $ (see ^Ts] above); for example, 8 = would be for radiation incident 
normal to the slab. 

21.6.1 Sample calculation in directory examples_exp/DSKRCTPBC 

Subdirectory examples_exp/DSKRCTPBC contains ddscat. par for calculating scattering by a 
0.0500/im thick Si3N4 slab supporting a doubly periodic array of Au disks, with periodicity 0.0800/im, 
disk diameter 0.0400/im, and disk height 0.0200/im, for light with wavelength A = 0.5320/xm. 

The TUC consists of a rectangular block of Si3N4, of dimension 15c? x 24(i x 24(i (8640 dipoles), 
supporting an Au disk of thickness 6d and a diameter ^12d (672 dipoles). With a "diameter" of only 
12d, the cross section of the "disk" is only roughly circular; each layer of the disk contains 1 12 dipoles. 

The volume of the ideal TUC is Vrvc = (0.08)^ x 0.05 + 7r(0.02)2 x 0.02 = 3.4513 x lO-^/^m^. 
Thus we set a^s = (3FTuc/47r)^/3 = 4.3514 x lO^^m- The DDA calculation has kd = .039377. 

The radiation is incident at an angle of 60° relative to the surface normal. The entire calculation 
required 2400 cpu sec on a 2.53 GHz cpu. Most of the cpu time was spent computing the effective A 
matrix for the calculation, which requires extensive summations; once this was obtained, the solution 
was found in 31 and 33 iterations, respectively, for the two incident polarizations, requiring ~8 cpu sec. 

21.7 HEXGONPBC = Target consisting of homogeneous hexagonal prism re- 
peated in target y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of hexagonal 
prisms. The individual prisms are assumed to be homogeneous and isotropic, just as for option RCTNGL 
(see |2aT8| ). 

Let us refer to a single hexagonal prism as the Target Unit Cell (TUC). The TUC is then repeated in 
the target y- and z-directions, with periodicities PYDxd and PZDxd, where d is the lattice spacing. To 
repeat in only one direction, set either PYD or PZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 

whereSHPARi, SHPAR2, SHPAR3, SHPAR4, SHPAR5 are numbers: 
SHPARi = prism length along prism axis (in units of d) in units of d 
SHPAR2 = 2xlength of one hexagonal side/d 

SHPAR3 = 1,2,3,4,5 or 6 to specify prism orientation in the TF (see below) 
SHPAR4 = PYD = periodicity in TF y direction/d 
SHPAR5 = PZD = periodicity in TF z direction/d 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
and SHPAR5. 

The physical size of the TUC is specified by the value of a^s (in physical units, e.g. cm), specified in 
the file ddscat .par, with the usual correspondence d — {in/SNy^^acs, where N is the number of 
dipoles in the TUC. 

With target option HEXGONPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions (assuming both NPY and NPZ are nonzero). 

The target axes (in the TF) are set to ai = xtf = (1iO,0)tf - i-c, normal to the "slab" - and 
a2 = Ytf = (0, 1,0)tf- 

The individual hexagons may have any of 6 different orientations relative to the slab: Let unit vectors 
h be II to the axis of the hexagonal prism, and let unit vector f be normal to one of the rectangular faces 



TARGET GENERATION: PERIODIC TARGETS 



of the hexagonal prism. Then 

SHPAR3=1 for h || Xtf, f |1 Ytf 
SHPAR3=2 for h || xtf, f || Ztf 
SHPAR3=3 for h II ytf, f || xtf 
SHPAR3=4 for h || yxF, f II Ztf 
SHPAR3=5 for h || Ztf, f || Xtf 
SHPAR3=6 for h || ztf, f |1 Ytf 



For example, one could construct a single infinite hexagonal column with the following line in 

ddscat . par: 



2.0 100.0 3 2.0 0. 



The TUC would be a thin hexagonal "slice" containing two layers of dipoles. The edges of the hexagon 
would be about 50d in extent, so the TUC would have approximately (3-\/3/2) x 50^ x 2 = 12990 
dipoles (13024 in the actual realization) within a 90 x 2 x 100 "extended target volume". The TUC 
would be oriented with the hexagonal axis in the yxF direction, with ztf normal to a rectangular faces 
of the prism (SHPAR3=3), and the structure would repeat in the yxF direction with a period of 2 x c? 
(SHPAR4=2 . 0). SHPAR5=0 means that there will be no repetition in the ztf direction. 

Note that SHPARi, SHPAR2, SHPAR4, and SHPAR5 need not be integers. However, SHPAR5, deter- 
mining the orientation of the prisms in the TF, can only take on the values 1,2,3,4,5,6. 

Important Note: For technical reasons, PYD and PZD must not be smaller than the "extended" 
target extent in the yxF and ztf directions. When the GPFAFT option is used for the 3-dimensional 
FFT calculations, the extended target volume always has dimensions/d = 2"3^5'^, where a, b, and c are 
nonnegative integers, with (dimension/d) > 1). 

The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles f3, Q, and $ (see 5 18 above); for example, 9 = would be for radiation incident 
normal to the periodic structure. 



21.8 LYRSLBPBC = Target consisting of layered slab, extended in target y and 
z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an array of multilayer bricks, layered in 
the Xtf dkection. The size of each brick in the yxF and ztf direction is specified. Up to 4 layers are 
allowed. 

The bricks are repeated in the yxF and Ztf direction with a specified periodicity. If Ly = Py and 
Lz = Pz, then the target consists of a continuous multilayer slab. For this case, it is most economical to 

set Ly/d = Lz/d = Py/d = Pz/d = 1. 

If Py = 0, then repetition in the yxF direction is suppressed - the target repeats only in the ztf direc- 
tion. 

If Pz — 0, then repetition in the Ztf direction is suppressed - the target repeats only in the yxF direction. 
The upper surface of the slab is asssume to be located at a;xF = 0. The lower surface of the slab is at 

Xtp = -Lx = - SHPARI Xrf. 

The multilayer slab geometry is specified with 9 parameters. The pertinent line in ddscat . par should 
read 



SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARe SHPAR7 SHPARg SHPARg 



where 

SHPARi = Lx/d = (brick thickness in xtf direction)/c? 
SHPAR2 = Ly/d = (brick extent in yTF direction)/^ 
SHPAR3 = Lz/d = (brick extent in zxf direction)/^ 
SHPAR4 = fraction of the slab with composition 1 
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SHPAR5 = fraction of the slab with composition 2 
SHPARg = fraction of the slab with composition 3 
SHPAR7 = fraction of the slab with composition 4 

SHPARg = P^j/d = (periodicity in ytf direction)/(i (0 to suppress repetition in yxF direction) 
SHPARg = Pz/d = (periodicity in ztf direction)/c? (0 to suppress repetition in ztf direction) 

For a slab with only one layer, set SHPAR5 = 0, SHPARg = 0, SHPAR7 = 0. 
For a slab with only two layers, set SHPARg = and SHPAR7 = 0. 
For a slab with only three layers, set SHPAR7 = 0. 

The user must set NCOMP equal to the number of nonzero thickness layers. 

The number N of dipoles in one TUC is iV = nint(SHPARi) x nint(SHPAR2) x nint(SHPAR3). 

The fractions SHPAR2, SHPAR3, SHPAR4, and SHPAR5 must sum to 1. The number of dipoles in 
each of the layers will be integers that are close to nint(SHPARi * SHPAR4), nint(SHPARi * SHPAR5), 
mnt(SHPARi * SHPARg), mnt(SHPARi * SHPAR7), 

The physical size of the TUC is specified by the value of a^s (in physical units, e.g. cm), specified 
in the file ddscat . par. Because of the way is defined (47raj^ff /3 = Nd'^), it should be set to 

aes = (3/47r)i/3ivV3rf= (3/47r)i/3(L,L,L,)i/3 



With target option LYRSLBPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i-C-, normal to the 
"slab" - and 3.2 = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles /?, Q, and $ (see ^Ts] above); for example, = would 
be for radiation incident normal to the slab. 

For this option, there are only two allowed scattering directions, corresponding to transmission and 
specular reflection. DDSCAT will calculate both the transmission and reflection coefficients. 
The last two lines in ddscat . par should appear as in the following example ddscat . par file. This 
example is for a slab with two layers: the slab is 26 dipole layers thick; the first layer comprises 76.92% 
of the thickness, the second layer 23.08% of the thickness. The wavelength is 0.532/im, the thickness is 

= (47r/3)i/3Ar2/3^^^, ^ (477/3)i/3(26)2/3o.o09189 = 0.1300Aim. 
The upper layer thickness is 0.2308L:r = 0.0300/im 

ddscat . par is set up to calculate a single orientation: in the Lab Frame, the target is rotated through 
an angle = 120°, with $ = 0. In this orientation, the incident radiation is propagating in the 
(—0.5, 0.866, 0) direction in the Target Frame, so that it is impinging on target layer 2 (Au). 

' =========== Parameter file for v7 . 2 ===================' 

PRELIMINARIES ****' 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 

'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 

'NOTBIN' = CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 

'**** Initial Memory Allocation 

26 1 1 = upper bounds on size of TUC 

'**** Target Geometry and Composition ****' 

'LYRSLBPBC = CSHAPE*9 shape directive 

26 1 1 0.7692 0.2308 11= shape parameters SHPARl - SHPAR9 

2 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/Eagle_2000' = refractive index 1 

' /u/draine/work/DDA/diel/Au_evap' = refractive index 2 

'**** Additional Nearfield calculation? ****' 

= NRFLD (=0 to skip nearfield calc., =1 to calculate nearfield E) 
(tract, extens . of calc. vol. in -x, +x, -y , +y , -z , +z ) 

Error Tolerance ****' 
l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC 1 E>-ACA | X>) / (NORM OF AC|E>) 
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'**** Maximum number of iterations ****' 
100 = MXITER 

'**** Interaction cutoff parameter for PBC calculations ****' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'**** Angular resolution for calculation of <cos>, etc. ****' 

0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] ) 

'**** Wavelengths (micron) ****' 

0.5320 0.5320 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 

Effective Radii (micron) **** ' 
0.009189 0.009189 1 'LIN' = eff. radii ( first , last , how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol. state eOl; =2 to also do orth. pol. state) 

'**** Specify which output files to write ****' 

I = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

120. 120. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

5 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 
'**** Specify Scattered Directions ****' 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
0. 0. = (M,N) for scattering 



21.9 RCTGL_PBC = Target consisting of homogeneous rectangular brick, ex- 
tended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of rectangular bricks. The 



bricks are assumed to be homogeneous and isotropic, just as for option RCTNGL (see 5 20. 18 1. 
Let us refer to a single rectangular brick as the Target Unit Cell (TUC). The TUC is then repeated 
in the j/tf- and ZTp-directions, with periodicities PYAEFFxoeff and PZAEFFXfloff, where Qcb = 
(3VTUc/47r)^/'^, where Vruc is the total volume of solid material in one TUC. This option requires 
5 shape parameters: 

The pertinent line in dd scat. par should read 



SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARi = (brick thickness)/(i in the xtf direction 
SHPAR2 = (brick thickness)/(i in the yTF direction 
SHPAR3 = (brick thickness)/(i in the ztf direction 
SHPAR4 = periodicity/d in the ?/tf direction 
SHPAR5 = periodicity/d in the ztf direction 



The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the j/tf and Ztf directions is determined by parameters 

SHPAR4 and SHPAR5. 

The physical size of the TUC is specified by the value of Ocff (in physical units, e.g. cm), specified in the 
file ddscat . par. 

With target option RCTGL_PBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i-C-, normal to the 
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"slab" - and 3.2 — (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles /?, Q, and $ (see ^Ts] above); for example, = would 
be for radiation incident normal to the slab. 



21.9.1 Sample calculation in directory examples_exp/RCTGL_PBC 

Subdirectory examples_exp/RCTGL_PBC contains ddscat . par for scattering by an infinite slab, 
constituted from 20 x 1 x 1 dipole TUCs. The wavelength A = 0.5/im, and the slab thickness is h ~ 
0.32/im. The slab has refractive index m = 1.50 + 0.02i. The interdipole spacing d — 0.32/im/20 = 
O.Oie/xm. The TUC has dimension y = 0.32/imx0.016^mx0.016Atm, andhence Ccff = (W/An)^^^ = 
(3 x 0.32 X 0.016 X 0.016/47r)i/3//m = 0.026942/im. The incident radiation is at an angle 9, = 40° 
relative to the surface normal. 



21.9.2 Sample calculation in directory examples_exp/RCTGL_PBC_NEARFIELD 



The directory examples_exp/RCTGL_PBC_NEARFIELD contains ddscat. par for calculation 
of scattering and absorption by an infinite slab of material with refractive index m = 1.50 + 0.02z and 
thickness h ~ 0.32/im in vacuo. The incident radiation has wavelength Avac = 0.5/j,m. The interdipole 



spacing is set to rf = h/20 = 0.016/im. This is an example problem shown in Fig. 7b of Draine & 
|Flatau| ( |2008] l. 

The slab is treated as a periodic array of Id x Id x 20d structures with periodicity Py = Id and 
= Id. The volume of the TUC is Vtuc = 20^^ 8.1920 x lO^^^m^, and the effective radius is 

fleff = (31/TUc/47r)i/3 = 0.026942/im. 



1.5 
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H.T, Draine 2013,01,24,1414 



Figure 9: Normalized electric field intensity |Ep/|Eop along a line normal to a slab of thickness h, refractive index m = 1.5 + 
0.02i. Radiation with X/h — 5is incident with incidence angle 9i = 40°. The interdipole spacing is set to d = h/20 — 0.005/im. 
This is the same problem reported in Fig. 7b of jDraine & Flatau|j2008) . 



21.10 RECRECPBC = Rectangular solid resting on top of another rectangular 
solid, repeated periodically in target y and z directions using periodic boundary 
conditions 

The TUC consists of a single rectangular "brick", of material 1, resting on top of a second rectangular 
brick, of material 2. The centroids of the two bricks along a line in the xtf direction. The bricks are as- 
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sumed to be homogeneous, and materials 1 and 2 are assumed to be isotropic. The TUC is then repeated 
in the ?/tf- and ztf -directions, with periodicities SHPAR4 x d and SHPAR5 x d. is the total volume of 
solid material in one TUC. This option requires 5 shape parameters; 
The pertinent line in dd scat. par should read 



SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 



Where 

SHPARi = (upper brick thickness)/(i in the xtf direction 
SHPAR2 = (upper brick thickness)/(i in the yxp direction 
SHPAR3 = (upper brick thickness)/(i in the ztf direction 
SHPAR4 = (lower brick thickness)/^ in the xtf direction 
SHPAR5 = (lower brick thickness)/ci in the yxF direction 
SHPARg = (lower brick thickness)/(i in the Ztf direction 
SHPAR7 = periodicity/d in the yxF direction 
SHPARg = periodicity/d in the ztf direction 



The actual numbers of dipoles Ni^, Niy, Niz, along each dimension of the upper brick, and N2X, 
^2y7 along each dimension of the lower brick, just be integers. Usually, N^j. — nint(SHPARi), 
Niy = nint(SHPAR2), A^iz = nint(SHPAR3), N2X = nint(SHPAR4), N2y = nint(SHPAR5), N2Z = 
nint(SHPAR6), where nint(a;) is the integer nearest to x, but under some circumstances Ni^, ^lyi -^iz, 
N2X, N2y, N2z might be larger or smaller by 1 unit. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi 
- SHPARg: 

N - (TVi, X Niy y.Niz) + (N2x x iVsy x (34) 

The periodicity in the yxF and Ztf directions is determined by parameters SHPAR7 and SHPARg. 
The physical size of the TUC is specified by the value of Ocff (in physical units, e.g. cm), specified in the 
file ddscat . par: 

d = {Air/mY'^a^s (35) 

The target is a periodic structure, of infinite extent in the yxF and ztf directions. The target axes are 
set to ai = Xtf - i e., normal to the "slab" - and 3.2 — yxF- The orientation of the incident radiation 



relative to the target is, as for all other targets, set by the usual orientation angles /3, 0, and <i> (see 5 18 
above) specifying the orientation of the target axes ai and 3.2 relative to the direction of incidence; for 
example, 9 = would be for radiation incident normal to the slab. 

The scattering directions are specified by specifying the diffraction order (Af , N); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix 5'-^'*-' calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is 

4 

(M, iV) =^4"'^/,„,j (36) 



where /i„ j is the incident Stokes vector See Draine & Flatau (20081 for interpretation of the Sij as 
transmission and reflection efficiencies. 



21.11 SLBHOLPBC = Target consisting of a periodic array of rectangular blocks, 
each containing a cylindrical hole 

Individual blocks have extent (a, fe, c) in the (xtf, yxF, ztf) directions, and the cylindrical hole has 
radius r. The period in the yxF-direction is Py, and the period in the zxF-direction is P^. 
The pertinent line in ddscat . par should consist of 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARg 
where SHPARi — a/ did is the interdipole spacing) 
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SHPAR2 = b/a 
SHPAR3 = c/a 
SHPAR4 = r/a 
SHPAR5 = Py/d 
SHPARg = Pz/d. 

Ideally, a/d =SHPARi, b/d =SHPAR2XSHPARi, and c/d ^SHPARaX SHPARi will be integers (so 
that the cubic lattice can accurately approximate the desired target). If Py — and P^ > the target is 
periodic in the z^f -direction only. 

If Pj^ > and Pz = the target is periodic in the yxF -direction only. 

If Pj, > it is required that Py > b, and if P^ > it is required that P^ > c, so that the blocks do not 
overlap. 

With Py = b and P^ = c, the blocks are juxtaposed to form a periodic array of cylindrical holes in a 
solid slab. 

Example: ddscat_SLBHOLPBC . par is a sample ddscat .par for a periodic array of cylindrical 
holes in a slab of thickness a, with holes of radius r, and period Py and Pd 

21.12 SPHRN_PBC = Target consisting of group of N spheres, extended in tar- 
get y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic array of A^— sphere structures, 
where one A^— sphere structure consists of N spheres, just as for target option NANSPH (see ^20.23[ ). 
Each sphere can be of arbitrary composition, and can be anisotropic if desired. Information for the de- 
scription of one A^-sphere structure is supplied via an external file, just as for target option NANSPH - 
see ^2023] ). 

Let us refer to a single A^-sphere structure as the Target Unit Cell (TUC). The TUC is then repeated 
in the i/tf- and zxp-directions, with periodicities PYAEFFxacfr and PZAEFFxacff, where Ooff = 
(3Vtuc/47I')^^'^, where Vruc is the total volume of solid material in one TUC. This option requires 
3 shape parameters: 

DIAMX = maximum extent of target in the target frame x direction/d 
PYAEFF = periodicity in target y direction/oefi 
PZAEFF = periodicity in target z direction/acg. 
The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 '^/enome' (quotes must be used) 

where 

SHPARi = target diameter in x direction (in Target Frame) in units of d 

SHPAR2= PYAEFF 
SHPAR3= PZAEFF. 

filename is the name of the file specifying the locations and relative sizes of the spheres. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameter SHPARi, 
which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. The 
physical size of the TUC is specified by the value of acs (in physical units, e.g. cm), specified in the file 
ddscat . par. 

The location of the spheres in the TUC, and their composition, is specified in file filename'. Please 
consult ^20. 23 [ above for detailed information concerning the information in this file, and its arrangement. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file filename ' should 
read 

a;i yi zi ri 1 1 1 
a;2 2/2 22 ra 1 1 1 
a:3 2/3 23 rs 1 1 1 



TARGET GENERATION: PERIODIC TARGETS 



i.e., every sphere has isotropic composition IC0MP=1 and the three dielectric function orientation angles 
are set to zero. 

When the user uses target option SPHRN_PBC, the target now becomes a periodic structure, of 
infinite extent in the target y- and z- directions. The target axis ai xlf ~ (1, 0, 0)tf - i-e., normal 
to the "slab" - and target axis a2 — ylf ~ (0, 1, 0)tf- The orientation of the incident radiation relative 
to the target is, as for aU other targets, set by the usual orientation angles /?, Q, and <i> (see 5 18 above). 



The scattering directions are, just as for other targets, determined by the scattering angles 9g, <j)g (see §22 
below). 

The scattering problem for this infinite structure, assumed to be illuminated by an incident monochro- 
matic plane wave, is essentially solved "exactly", in the sense that the electric polarization of each of the 
constituent dipoles is due to the electric field produced by the incident plane wave plus all of the other 
dipoles in the infinite target. 

However, the assumed target will, of course, act as a perfect diffraction grating if the scattered radia- 
tion is calculated as the coherent sum of all the oscillating dipoles in this periodic structure: the far-field 
scattered intensity would be zero in all directions except those where the Bragg scattering condition is 
satisfied, and in those directions the far-field scattering intensity would be infinite. 

To suppress this singular behavior, we calculate the far-field scattered intensity as though the sepa- 
rate TUCs scatter incoherently. The scattering efficiency Qsca and the absorption efficiency Qabs are 
defined to be the scattering and absorption cross section per TUC, divided by Tra^g, where Ocff = 
(3VTUc/47r)^/'^, where Vruc is the volume of solid material per TUC. 

Note: the user is allowed to set the target periodicity in the target i/tf (or ztf) direction to values that 
could be smaller than the total extent of one TUC in the target ?/tf (or ztf) direction. This is physically 
allowable, provided that the spheres from one TUC do not overlap with the spheres from neighboring 
TUCs. Note that DDSCAT does not check for such overlap. 



21.12.1 Sample calculation in directory examples_exp/SPHRN_PBC 

Subdirectory example s_exp/SPHRN_PBC contains ddscat . par to calculate scattering by a doubly- 
periodic array with the target unit cell consisting of a random cluster of 16 spheres. The calculation is 
carried out with double precision arithmetic. Because convergence is slow, the error tolerance is set 
to TOL = 5 . e-5 rather than the usual 1 . e-5, and the maximum number of iterations allowed is 
increased to MX ITER = 2 000. 



21.13 TRILYRPBC = Three stacked rectangular blocks, repeated periodically 

The target unit cell (TUC) consists of a stack of 3 rectangular blocks with centers on the xtf axis. The 
TUC is repeated in either the ytf direction, the ztf direction, or both. A total of 11 shape parameters 
must be specified: 

SHPARi = x-thickness of upper layer/d [material 1] 

SHPAR2 = y-width/d of upper layer 

SHPAR3 = z-width/d of upper layer 

SHPAR4 = x-thickness/d of middle layer [material 2] 

SHPAR5 = y-width/d of middle layer 

SHPARg = z-width/d of middle layer 

SHPAR7 = x-width/d of lower layer 

SHPARg = y-width/d of lower layer 

SHPARg = z-width/d of lower layer 

SHPARio = period/d in y direction 

SHPARii= period/d in z direction 



SCATTERING DIRECTIONS 



22 Scattering Directions 
22.1 Isolated Finite Targets 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported 
in the output files wxxxryyykzzz . sea . The scattering direction is specified through angles Og and ipg 
(not to be confused with the angles 8 and $ which specify the orientation of the target relative to the 
incident radiation!). 

For isolated finite targets (i.e., PBC not employed) there are two options for specifying the scatter- 
ing direction, with the option determined by the value of the string CMDFRM read from the input file 
ddscat . par. 

1. If the user specifies CMDFRM=' LFRAME' , then the angles 0, (j) input from ddscat .par are 
understood to specify the scattering directions relative to the Lab Frame (the frame where the 
incident beam is in the x— direction). 

When CMDFRM=' LFRAME' , the angle is simply the scattering angle Og'. the angle between the 
incident beam (in direction xlf) and the scattered beam {Og = for forward scattering, Og — 180° 
for backscattering). 

The angle specifies the orientation of the "scattering plane" relative to the xlf — Ylf plane. 
When = the scattering plane is assumed to coincide with the xlf — Ylf plane. When — 90° 
the scattering plane is assumed to coincide with the xlf — zlf plane. Within the scattering plane 
the scattering directions are specified by < 6* < 180°. Thus: 

As = Xlf cos 6 + ylf sin 9 cos 4> + zlf sin 6 sin 4> , (37) 



2. If the user specifies CMDFRM=' TFRAME' , then the angles 9, 4> input from ddscat .par are 
understood to specify the scattering directions fis relative to the Target Frame (the frame defined 
by target axes xtf, Ytf, ztf)- S is the angle between and xtf, and (p is the angle between the 
hg — xtf plane and the xtf — Ytf plane. Thus: 

= xxf cos + yxF sin cos (/) + zxf sin sin (/) . (38) 



Scattering directions for which the scattering properties are to be calculated are set in the parameter 
file ddscat . par by specifying one or more scattering planes (determined by the value of (pg) and for 
each scattering plane, the number and range of 9s values. The only limitation is that the number of 
scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code as distributed it is set 
to MXSCA=1000). 



22.2 Scattering Directions for Targets that are Periodic in 1 Dimension 



For targets that are periodic, scattering is only allowed in certain directions ( Draine & Flatau|2008 see). 



If the user has chosen a PBC target (e.g, CYLNDRPBC, HEXGONPBC, or RCTGL_PBC), SPHRN_PBC 
but has set one of the periodicities to zero, then the target is periodic in only one dimension - e.g., 
CYLNDRPBC could be used to construct a single infinite cylinder 

In this case, the scattering directions are specified by giving an integral diffraction order AI = 
0,±1,±2,... and one angle, the azimuthal angle around the target repetition axis. The diffraction 
order M determines the projection of onto the repetition direction. For a given order M, the scatter- 
ing angles with = — > 27r form a cone around the repetition direction. 

For example, if PYD > (target repeating in the ?/tf direction), then M determines the value of 
kgy = ks • yTF, where yxF is the unit vector in the Target Frame direction: 



2nM/Ly (39) 



where k^y = kg • yxF. where kg is the incident k vector Note that the diffraction order M must satisfy 
the condition 

[koy - ko){Ly/2T:) < M < (fco - koy){Ly/2T:) (40) 



SCATTERING DIRECTIONS 



where Ly = PYD x d is the periodicity along the y axis in the Target Frame. M = is always an allowed 
diffraction order. 

The azimuthal angle C. defines a right-handed rotation of the scattering direction around the target 
repetition axis. Thus for a target repetition axis yxF. 

ksx = fc_LC0sC , (41) 
ksz = k±smC , (42) 

koy + 2'jTM/Ly. For a target with repetition axis ^tf. 

ksx = fc_LC0sC , (43) 
ksy = fc_LsinC , (44) 

where kx = {kp - fc^J/Z^, ksz = koz + 2nM/Lz. 

The user selects a diffraction order M and the azimuthal angles ( to be used for that M via one line 
in ddscat . par. An example would be to use CYLNDRPBC to construct an infinite cyUnder with the 
cylinder direction in the yxF direction: 

' ============== Parameter file for v7.2 ===================' 

Preliminaries 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, GPBICG, QMRCCG, PETRKP) — solution method 
'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) — FFT implementation 
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 

Initial Memory Allocation 
100 100 100 = size allowance for target generation 

Target Geometry and Composition 
'CYLNDRPBC = CSHAPE*9 shape directive 

1 64.499 2 1.0 0.0 = shape parameters SHPARl , SHPAR2, SHPAR3, SHPAR4, SHPAR5, 

1 = NCOMP ^ number of dielectric materials 

' . . /diel/ml . 33_0 . 1 ' ^ name of file containing dielectric function 
Additional Nearfield calculation? 

= NRFLD (=0 to skip nearfield calc, =1 to calculate nearfield E) 

0.0 0.0 0.0 0.0 0.0 0.0 (tract, extens. of calc. vol. in -x, +x, -y, +y, -z, +z) 
' **** Error Tolerance 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

Maximum number of iterations 
100 = MXITER 

Interaction cutoff parameter for PBC calculations 
5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

Angular resolution for calculation of <cos>, etc. ****' 
0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] ''2 ) 
' -k-k-k-k Vacuum wavelengths (micron) 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 

' ±kk-k Refractive index of ambient medium' 
1.0000 = NAMBIENT 

'★*** Effective Radii (micron) **** ' 

2.8555 2.8555 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol. state) 
'k-k-k-k Specify which output files to write ****' 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient.) 
'**** Prescribe Target Rotations 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
0. 0. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
' *kkk Specify first IWAV, IRAD, lORI (normally 0) 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
' kkkk Specify Scattered Directions ****' 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering cones 

0. 0. 180. 0.05 = M, zeta_min, zeta_max, dzeta (in degrees) 



where k± = {k^ - kly)^/"^, with ksy = 



SCATTERING DIRECTIONS 



In this example, a single diffraction order M = is selected, and is to run from Cmin = to 
Cmax = 180° in increments of = 0.05°. 

There may be additional lines, one per diffraction order. Remember, however, that every diffraction 
order must satisfy eq. (40 1. 



22.3 Scattering Directions for Targets for Doubly-Periodic Targets 

If the user has specified nonzero periodicity in both the y and z directions, then the scattering direc- 
tions are specified by two integers - the diffraction orders M, N for the y, z directions. The scattering 
directions are 



ko 



2ttM 



ks = ± — (ko • xtf) xtf + ko • yxF + — — ) Ytf + ( ko • ztf + — — ) ztf 



J 



2nN 



J 



(45) 



fco - ko • ytf 



27rMy 



ko • Ztf 



2ttNY 



-I 1/2 



(46) 



where the + sign gives transmission, and the — sign gives reflection. The integers M and N must 
together satisfy the inequality 



(ko • YTF + 2TiM/Lyf + (ko • ztf + '^^N/L.f < 



(47) 



which, for small values of Ly and L^, may limit the scattering to only (M, N) — (0, 0). [Of course, 
(0,0) is always allowed]. For each {M,N) specifi ed in ddscat .par, DDSCAT 7.2 will calculate 
the generalized Mueller matrix 5'|J''' (M, A^) - see 5 26.2 At large distances from the infinite slab, the 
scattered Stokes vector in the (M, TV) diffraction order is 



LcaAM, N) 



4 

E 



SlfiM,N)hr.,j 



(48) 



where lin.j is the incident Stokes vector, and S^j'^^ (M, N) is the generalization of the 4x4 Miiller 
scattering matrix to targets that are periodic in 2-directions ( Draine & Flatau|2008| ). There are distinct 



(M, N) for transmission and for reflection, corresponding to the ± in eq. (45 i 



Draine & Flatau 
mission coethcients 



(2008|l (eq. 69-71) show how the sf^'^\M, N) are easily related to familiar "trans- 



and "reflection coefficients". 
Here is a sample ddscat . par file as an example: 



' NOTORQ' 
■ PBCGS2' 
' GPFAFT' 
' GKDLDR' 
■ NOTBIN' 



Parameter file for v7.2 
Preliminaries 

(DOTORQ, 
(PBCGS2, 
(GPFAFT, 
(GKDLDR, 
(NOTBIN, 



CMT0RQ*6 
CMDS0L*6 
CMETHD*6 
CALPHA*6 
CBINFLAG 



— either do or skip torque calculations 
GPBICG, QMRCCG, PETRKP) — solution method 

— FFT implementation 

— recipe for dipole polarizabilities 
ALLBIN) — binary output? 



1 



NOTORQ 
PBCGST 
FFTMKL 
LATTDR 
ORIBIN 

**** Initial Memory Allocation 

00 100 100 - dimensioning allowance for target generation 
**** Target Geometry and Composition 
RCTGL_PBC' = CSHAPE*9 shape dirctive 

0.0 1.0 1.0 1.0 1.0 = shape parameters PARI, PAR2, PAR3,... 

^ NCOMP ^ number of dielectric materials 
. . /diel/ml . 33_0 . 1 ' = name of file containing dielectric function 

Additional Nearfield calculation? **** 
= NRFLD (^0 to skip nearfield calculation, -1 to calculate nearfield E) 
.0 0.0 0.0 0.0 0.0 0.0 (tract, extens . of calc. vol. in -x, +x, -y, +y, -z, +z 
**** Error Tolerance v^^^v^v^' 

.OOe-5 = TOL = max allowed (norm of | G>=AC I E>-ACA | X>) / (norm of AC|E>) 
**** Maximum number of iterations 
00 = MXITER 

■k-k-k-k Integration limiter for PBC calculations 



INCIDENT POLARIZATION STATE 



5.00e-3 ^ GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

Angular resolution for calculation of <cos>, etc. 
0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] "2 ) 

Wavelengths (micron) 
6.283185 6.283185 1 ' INV = wavelengths ( fir st , last , how many , how=LIN, INV, LOG) 

Effective Radii (micron) ' 
2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
'★*** Define Incident Polarizations ***★' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pel. state eOl; =2 to also do orth . pol . state) 

'**** Specify which output files to write ****' 

I = IWRKSC (-0 to suppress, ^1 to write ".sea" file for each target orient.) 
' **** Specify Target Rotations *+++' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

60. 60. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 - indices ij of elements to print 
'**** Specify Scattered Directions ****' 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
= M, N (diffraction orders) 



23 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the xlf 
axis. DDSCAT allows the user to specify a general elliptical polarization for the incident radiation, by 
specifying the (complex) polarization vector egi. The orthonormal polarization state eo2 = xlf x 
is generated automatically if ddscat . par specifies I0RTH=2. 

For incident linear polarization, one can simply set egi = y by specifying 
(0,0) (1,0) (0,0) 

in ddscat .par; then §02 = z- For polarization mode eoi to correspond to right-handed circular po- 
larization, set Gqi = (y + iz)/-y2 by specifying (0,0) (1,0) (0,1) in ddscat . par (DDSCAT 
automatically takes care of the normalization of Gqi); then §02 — {iy + z)/y/2, corresponding to left- 
handed circular polarization. 



24 Averaging over Scattering Directions: g{l) = (cos 6s), etc. 
24.1 Angular Averaging 



An example of scattering by a nonspherical target is shown in Fig. 10 showing the scattering for a tilted 
cube. Results are shown for four different scattering planes. 

DDSCAT automatically carries out numerical integration of various scattering properties, including 



(cos 0s 

.2 a 



• (cos' 

• g = (cos 6's)xlf + (sin 9s cos < 



)y + (sin^s sin(/)s)z (see516i; 



Qr, provided option DOTORQ is specified (see 5 16 1 



The angular averages are accomplished by evaluating the scattered intensity for selected scattering di- 
rections {9s, ips), and taking the appropriately weighted sum. Suppose that we have Ng different values 
of 9s, 

9 = 9j, ]^l,...,Ne , (49) 
and for each value of 9j, N^{j) different values of (f>s: 



(t>s^(t>jM, k = l,...,N^{j) 



(50) 
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30 60 90 120 150 180 

es(deg) B,I,D„,„=3003 „,0; 



Fig lire 10; Scattered intensity for incident unpolarized light for a cube witii m = 1.5 + O.OOli and x = 2naefs/X ~ 6.5 (i.e., 
d/X = 1.6676, where d is the length of a side). The cube is tilted with respect to the incident radiation, with O = 30°, and rotated 
by /3 — 15° around its axis to break reflection symmetry. The Mueller matrix element Sn is shown for 4 different scattering planes. 
The strong forward scattering lobe is evident. It is also seen that the scattered intensity is a strong function of scattering angle (ps as 
well as ^3 ~ at a given value of 8s (e.g., 9^ = 120°), the scattered intensity can vary by orders of magnitude as (jts changes by 90°. 



For a given j, the values of (pj^k are assumed to be uniformly spaced: (j)j^k+i ^ (pj.k = 27r/A^0(j). The 
angular average of a quantity f{6s, (ps) is approximated by 



2ir 



d4)sf{0s,(t>s 



47r 



j=i fc=i 



Ne,k 



[cos(0j_i) - cos(0j+i)] , j = 2,...,iVe-l 



27r 
2tt 



For a sufficiently large number of scattering directions 

Ng 



cos(6'i) + cos(6'2) 



jONg-l) +cos{9Ng) ^ 
2 



^,ea = E^^(^') 



(51) 

(52) 
(53) 
(54) 

(55) 



the sum (51 1 approaches the desired integral, but the calculations can be a significant cpu-time burden, 



so efficiency is an important consideration. 



24.2 Selection of Scattering Angles 9s, (ps 

Beginning with DDSCAT 6.1, an improved approach is taken to evaluation of the angular average. Since 
targets with large values of x = 27raoff /A in general have strong forward scattering, it is important to 
obtain good sampling of the forward scattering direction. To implement a preferential sampling of the 
forward scattering directions, the scattering directions {9, (p) are chosen so that 9 values correspond 
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to equal intervals in a monotonically increasing function s{9). The function s{d) is chosen to have a 
negative second derivative (Ps/d9^ < so that the density of scattering directions will be higher for 
small values of 9. DDSCAT takes 



(56) 



This provides increased resolution (i.e., increased ds/dO) for 9 < 9q. We want this for the forward 
scattering lobe, so we take 

^0 = . (57) 

1 + X 

The 9 values run from 6 — to 9 — tt, corresponding to uniform increments 

1 



As = I TT 

Ng-1\ 7T + Oo 



If we now require that 



c[A6i] 



As 7r/2 
{ds/d9)g^T^ 3 + a; 



we determine the number of values of 9: 



Ng = l + 



2(3 + x) [l + l/(7r + go)] 
11 [l + 9^/{n + 9oY] 



(58) 



(59) 



(60) 



Thus for small values of x, max[A6'] — 30° t], and for x ^ 1, max[A6'] — > %{)°rj/x. For a sphere, 
minima in the scattering pattern are separated by ^ 180° /a;, so 77 = 1 would be expected to marginally 
resolve structure in the scattering function. Smaller values of -q will obviously lead to improved sampling 
of the scattering function, and more accurate angular averages. 
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sphere: m= 1.33 + 0. Oli, N=17904 

, \ , , , \ , , ^ 

= 2 |error[(cos6)] 




fi T Draine 2003.1 1 0' 



2 x = 27Ta/X 4 

Figure 1 1 : Errors in (cos 9) and (cos^ 9) calculated for a A'^ = 17904 dipole pseudosphere with m = 1.33 + O.Oli, as functions 
of a:: = 27ra/A. Results are shown for different values of the parameter 77. 



The scattering angles 9j used for the angular averaging are then given by 



{s,-l-9a)+[{l + 9a-s,f+A9at 



,1/2 



(61) 
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Figure 12: Same as Fig.[TT] but for a iV = 32768 dipole cube with m = 1.5 + O.OOli. 



where 



1 



3 



l,...,iVe 



(62) 



For each Bj, we must choose values of 0. For 6*1 = and O^e ^ only a single value of (p is needed 
(the scattering is independent of </) in these two directions). For < 6j < tt we use 



= max {3, nint [47r sin(6'j)/(6'j+i - 6'j-i)]} , 
where nint = nearest integer. This provides sampling in cj) consistent with the sampling in 9. 



(63) 



24.3 

Figure 



Accuracy of Angular Averaging as a Function of rj 

shows the absolute errors in (cos 0) and (cos^ &) calculated for a sphere with refractive index 
O.Oli using the above prescription for choosing scattering angles. The error is shown as a 

1, and that 



11 



m = 1.33 

function of scattering parameter x. We see that accuracies of order 0.01 are attained with -q 
the above prescription provides an accuracy which is approximately independent of x. We recommend 
using values of 77 < 1 unless accuracy in the angular averages is not important. 



25 Scattering by Finite Targets: The Mueller Matrix 
25.1 Two Orthogonal Incident Polarizations (1 0RTH=2) 

Throughout the following discussion, x, y, z are unit vector defining the Lab Frame (thus ko — x). 
DDSCAT internally computes the scattering properties of the dipole array in terms of a complex 



scattering matrix fmi{9s, (ps) (Draine 1988) , where index I ~ 1,2 denotes the incident polarization state. 



m = 1, 2 denotes the scattered polarization state, and 9s,(t>s specify the scattering direction. Normally 
DDSCAT is used with I0RTH=2 in ddscat .par, so that the scattering problem will be solved for 
both incident polarization states (1^1 and 2); in this subsection it will be assumed that this is the case. 
Incident polarization states / = 1, 2 correspond to polarization states Gqi, 602; recall that polarization 
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state eoi is user-specified, and eo2 = x x Gq^ " Scattered polarization state m = 1 corresponds to 
linear polarization of the scattered wave parallel to the scattering plane (ei = e||s = dg) and m = 2 
corresponds to linear polarization perpendicular to the scattering plane (in the +0^ direction). The 
scattering matrix /,,„; was defined Draine (J988 ) so that the scattered electric field is related to the 
incident electric field (0) at the origin (where the target is assumed to be located) by 



Es- 
E. • 



cxp(zk, ■ r) / /ii 

kr V /21 /22 



E.(0)-e5i 



(64) 



The 2x2 complex scattering amplitude matrix (with elements 5*1, S2, S3, and S4) is defined so that 
( IBohren & Huffman] 1983| see) 



E« 
-E, 



expjiks ■ r) / 5*2 5*3 
—ikr \ S4 Si 



E,(0)-e,|| 
E.(0) •e.i 



(65) 



where e^n , e^j^ are (real) unit vectors for incident polarization parallel and perpendicular to the scattering 
plane (with the customary definition of e^j^ — x xlf)- 



From ( 64|65 1 we may write 



'5'2 S3 

S4 Si 



E,(0)-e,|| 
E.(0) •e.i 



hi 

-hi 



/12 

— f22 



E.(0)-e5i 
E,(0)-e52 



(66) 



Let 



a — 


^01 




(67) 


b = 


^01 


zlf , 


(68) 


c = 


^02 


YLF , 


(69) 


d = 


^02 


zlf • 


(70) 



Note that since eoi, eo2 could be complex (i.e., elliptical polarization), the quantities a, &, c, d are com- 
plex. Then 



-01 
-02 



a b 
c d 



YLF 

Zlf 



and eq. ( 66 1 can be written 
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S4 Si 
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E.(0) •e.i 



— fll —fl2 
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The incident polarization states eq and e^j^ are related to xlf, Ylf, zlf by 



YLF 

Zlf 



cos ( 
sin ( 



sm (p. 

- COS( 



e.,|| 



(71) 



(72) 



(73) 



substituting (73 1 into ( [72| we obtain 

52 S3\f E,(0) • e,j| \ Y -fii -fi2 \f a 6 \ / cos0, sin0, \ / E,(0) • e,,, 
54 Si )\ E,(0) -e,^ ) '\ hi h2 )\c d )\ sin0, -cos0, ) \ E,(0) • e,j_ 

(74) 

Eq. (74 1 must be true for all Ei(0), so we obtain an expression for the complex scattering amplitude 
matrix in terms of the fmi '■ 



S2 S3 
S4 Si 



— fll —fl2 
hi /22 



a b 
c d 



cos ( 
sin f 



sm (p. 
~ cos ( 



(75) 



Draine 



1 1988 1 adopted the convention eo2 = x x e*)i. The customary definition of is e^^ = e^ii x k. Thus, if eoi — 



then eo2 = — ei±. 

A frequent choice is eoi = Ylp, with eo2 ~ xlf x ylf = zlf- 
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This provides the 4 equations used in subroutine GETMUELLER to compute the scattering amplitude 
matrix elements: 

51 = -i [/2i(6cos(/)s - asin(/)s) + /22(dcos(/)s - csin(/)s)] , (76) 

52 = -i[fii{acos(t)s + bsm(t)s) + fi2{ccos(l)s+dam(t)s)] , (77) 

53 = i[fii{bcos(f>s - asm(j>s) + fi2idcos(l>s - csm(l>s)] , (78) 
84 = i [/2i(acos0s + 6sin0s) + /22(ccos0s + dsin^s)] . (79) 



25.2 Stokes Parameters 

It is both convenient and customary to characterize both incident and scattered radiation by 4 "Stokes 
parameters" - the elements of the "Stokes vector". There are different conventions in the literature; we 
adhere to the definitions of the Stokes vector (I,Q,U,V) adopted in the excellent treatise by ,Bohren & 
Huffman] ( |1983| l, to which the reader is referred for further detail. Here are some examples of Stokes 
vectors (J, Q, U, V) = (1, Q//, U/I, V/I)I: 

• (1, 0, 0, 0)/ : unpolarized light (with intensity /); 

• (1, 1, 0, 0)/ : 100% linearly polarized with E parallel to the scattering plane; 

• (1, —1, 0, 0)/ : 100% linearly polarized with E perpendicular to the scattering plane; 

• (1, 0, 1, 0)/ : 100% linearly polarized with E at +45° relative to the scattering plane; 

• (1, 0, —1, 0)/ : 100% linearly polarized with E at -45° relative to the scattering plane; 

• (1,0, 0, 1)/ : 100% right circular polarization {i.e., negative helicity); 

• (1,0, 0, —1)/ : 100% left circular polarization (i.e., positive helicity). 



25.3 Relation Between Stokes Parameters of Incident and Scattered Radiation: 
The Mueller Matrix 

It is convenient to describe the scattering properties of a finite target in terms of the 4x4 Mueller ma- 
trix Sij relating the Stokes parameters (J^, Q^, Ui, Vi) and {Is^ Qs, Us, Vs) of the incident and scattered 
radiation: 
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(80) 
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Once the amplitude scattering matrix elements are obtained, the Mueller matrix elements can be com- 
puted ( |Bohren & Huffman|1983| l: 

\S4\')/2 , 
\S3n/2 , 
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(81) 

These matrix elements are computed in DDSCAT and passed to subroutine WRITE SCA which handles 
output of scattering properties. Although the Muller matrix has 16 elements, only 9 are independent. 

The user can select up to 9 distinct Muller matrix elements to be printed out in the output files 
wxxxryyykzzz .sea and wxxxryyyori . avg; this choice is made by providing a list of indices in ddscat 
(see Appendix [a). 

If the user does not provide a list of elements, WRITESCA will provide a "default" set of 6 selected 
elements: ^n, ^21, 5*31, >5'4i (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), 512, and ^la. 

In addition, WRITESCA writes out the linear polarization P of the scattered light for incident unpo- 
larized light ( Bohren & Huftman|1983| see): 



P = 



(Si 



8i\ 



(82) 



par 



25.4 Polarization Properties of the Scattered Radiation 

The scattered radiation is fully characterized by its Stokes vector (/s, Qs, f^s, ^s)- As discussed in 
[Bohren & Huffman] ( |1983p (eq. 2.87), one can determine the linear polarization of the Stokes vector 
by operating on it by the Mueller matrix of an ideal linear polarizer: 
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(83) 



where ^ is the angle between the unit vector 0^ parallel to the scattering plane ("SP") and the "transmis- 
sion" axis of the linear polarizer. Therefore the intensity of light polarized parallel to the scattering plane 
is obtained by taking ^ = and operating on (J^, Qs,Us:Vs) to obtain 



/(E, II SP) 



1 



-Sis 



2kh 



[(^U + 82l)h + {8,2 + 822)Q^ + {813 + S23)U, + 



(84) 

(5i4 + ^24)^,] . (85) 
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Similarly, the intensity of light polarized perpendicular to the scattering plane is obtained by taking 

€ = 7r/2: 

/(E,±SP) = \{Is-Qs) (86) 

= - S2l)h + (5i2 - S22)Qi + (5i3 - S23)U, + (5i4 " S^^Wi] . (87) 

25.5 Relation Between Mueller Matrix and Scattering Cross Sections 

Differential scattering cross sections can be obtained directly from the Mueller matrix elements by noting 
that 



1 / dCsca. 

dQ. 



li . (88) 

Let SP be the "scattering plane": the plane containing the incident and scattered directions of propaga- 
tion. Here we consider some special cases: 

• Incident Ught unpolarized: Stokes vector Si = /(1, 0, 0, 0): 

- cross section for scattering with polarization Eg || SP: 

'^■^(IS.i^ + lftl'j-ilSn+S.) (89, 

- cross section for scattering with polarization _L SP: 

- total intensity of scattered Ught: 

• Incident light polarized with Ej || SP: Stokes vector Sj = 7(1, 1, 0, 0): 

- cross section for scattering with polarization Eg || to SP: 

- cross section for scattering with polarization Eg ± to SP: 

^ = ^l^^l' = i (^11 + - - S22) (93) 

- total scattering cross section: 

^ = ^(|52p + 1^4^ = ^(511 + 5,.) (94) 

• Incident Ught polarized with Ej _L SP: Stokes vector Si = 7(1, —1, 0, 0): 

- cross section for scattering with polarization Eg || to SP: 

- cross section for scattering with polarization E^ _L SP: 
— = ^l^'^il' = ^2 (^11 - - S21 + S22) (96) 
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- total scattering cross section: 



"'^^'^ 4 {\Si\' + IS^n = ^ (^11 - S,,) (97) 



• Incident light linearly polarized at angle 7 to SP: Stokes vector = /(I, cos 27, sin 27, 0) 
- cross section for scattering with polarization || to SP: 



df! " 2fc2 t^*^" + ^21) + (5i2 + ^22) cos 27 + (5i3 + ^23) sin 27] (98) 

- cross section for scattering with polarization _L SP: 

dC 1 

= ^ [{Sii -821) + (5i2 - 522) cos 27 + (5i3 - ^23) sin 27] (99) 

- total scattering cross section: 

^ = ^ [Sn + S,2 cos 27 + ^13 sin 27] (100) 

25.6 One Incident Polarization State Only (I0RTH=1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state - for 
example, when each solution is very time-consuming, and the target is known to have some symmetry 
so that solving for a single incident polarization state may be sufficient for the required purpose. In this 
case, set I0RTH=1 in ddscat . par. 

When I0RTH=1, only /n and /21 are available; hence, DDSCAT cannot automatically generate 
the Mueller matrix elements. In this case, the output routine WRITE SCA writes out the quantities 
1/21 p, Re(/ii/|]^), and Im(/ii/|]^) for each of the scattering directions. 

The differential scattering cross section for scattering with polarization Eg || and _L to the scattering 
plane are 

'^^'^'^ ^ 1/11 P (101) 



dn fc2' 

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 even if 
ddscat . par specified I0RTH=1: this is because when more than one value of the target orientation 
angle $ is required, there is no additional "cost" to solve the scattering problem for the second incident 
polarization state, since when solutions are available for two orthogonal states for some particular target 
orientation, the solution may be obtained for another target orientation differing only in the value of 
$ by appropriate linear combinations of these solutions. Hence we may as well solve the "complete" 
scattering problem so that we can compute the complete Mueller matrix. 



26 Scattering by Periodic Targets: Generalized Mueller Matrix 



The Mueller scattering matrix Sij{6) described above was originally defined (Bohren & Huffman 1983 
see, e.g.,) to describe scattering of incident plane waves by finite targets, such as aerosol particles or dust 
grains. 



Draine & Flatau ( 2008|l have extended the Mueller scattering matrix formalism to also apply to 



targets that are periodic and infinite in one or two dimensions (e.g., an infinite chain of particles, or a 
two-dimensional array of particles). 
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26.1 Mueller Matrix S'f^'^^ for Targets Periodic in One Direction 



The dimensionless 4x4 matrix S^^'^'' (Af, () describes the scattering properties of targets that are periodic 
in one dimension. In the radiation zone, for incident Stokes vector /jn, the scattered Stokes vector in 
direction (M, C) (see ^22.2[ ) is (see |Draine & FlataupOOS] eq. 63) 



(103) 



where R is the distance from the target repetition axis. 



DDSCAT 7.2 reports the scattering matrix elements S'^"''-' in the output files waaarbbbkccc .sea 



26.2 Mueller Matrix S^]"^^ (M, A^) for Targets Periodic in Two Directions 

For targets that are periodic in two directions, the scattering intensities are described by the dimensionless 
4x4 matrix sI'^'^\m, N). In the radiation zone, for incident Stokes vector the scattered Stokes 



vector in scattering order (M, N) is (see Draine & Flatau 2008 eq. 64) 



/,ca,^(M, iV) = V Slf{M, N) Jin,,- 



(104) 



where integers (A/, N) define the scattering order (see §22.3 1. For given {A-I, N) there are actually two 
possible scattering directions - corresponding to transmission and reflection (see ^22.3[ l, and for each 
there is a value of S'fJ'*^ (M, N), which we will denote S'fJ''' (M, N, tran) and S'||'^^ (M, N, refl). 

For targets with 2-d periodicity, the scattering matrix eleme nts iS'^J'^'' (M, N) are di rectly related to 
the usual transmission coefficient and reflection coefficient (see [Draine & Flatau|2008| eq. 69-71). For 
unpolarized incident radiation, the reflection coefficient R, transmission coefficient T, and absorption 
coefficient A are just 



M,N 

R = ^5(f)(A^,Ar,refl) 

M,N 

A = l-T-R . 



(105) 

(106) 
(107) 



27 Composite Targets with Anisotropic Constituents 

Section|20]includes targets composed of anisotropic materials (ANIELLIPS, ANIRCTNGL, ANI_ELL_2, 
ANI_ELL_3). However, in each of these cases it is assumed that the dielectric tensor of the target ma- 
terial is diagonal in the "Target Frame". For targets consisting of a single material (in a single domain), 
it is obviously possible to choose the "Target Frame" to coincide with a frame in which the dielectric 
tensor is diagonal. 

However, for inhomogeneous targets containing anisotropic materials as, e.g., inclusions, the optical 
axes of the constituent material may be oriented differently in different parts of the target. In this case, it 
is obviously not possible to choose a single reference frame such that the dielectric tensor is diagonalized 
for all of the target material. 

To extend DDSCAT to cover this case, we must allow for off-diagonal elements of the dielectric 
tensor, and therefore of the dipole polarizabilities. It will be assumed that the dielectric tensor e is 
symmetric (this excludes magnetooptical materials - check this). 

Let the material at a given location in the grain have a dielectric tensor with complex eigenvalues e^^\ 
j — 1 — 3. These are the diagonal elements of the dielectric tensor in a frame where it is diagonalized 
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(i.e., the frame coinciding with the principal axes of the dielectric tensor). Let this frame in which the 
dielectric tensor is diagonalized have unit vectors corresponding to the principal axes of the dielectric 
tensor We need to describe the orientation of the "Dielectric Frame" (DF) - defined by the "dielectric 
axes" Bj - relative to the Target Frame. It is convenient to do so with rotation angles analogous to the 
rotation angles Q, ^, and (3 used to describe the orientation of the Target Frame in the Lab Frame: 
Suppose that we start with unit vectors aligned with the target frame Xj . 

1. Rotate the DF through an angle 6'df around axis yxF, so that 6'df is now the angle between ei 
and XTp. 

2. Now rotate the DF through an angle (/)df around axis xtf, in such a way that 62 remains in the 
xtf — Gi plane. 

3. Finally, rotate the DF through an angle /3df around axis ei. 

The unit vectors are related to the TF basis vectors Xtf, Ytf. ^tf by" 

ei = Xtf cos 6'df + yxF sin 6'df cos (/>DF + ztf sin 6'df sin (/)DF (108) 

62 = —Xtf sin ^df cos /3df + Ytf [cos 6'df cos /3df cos (/)df — sin /3df sin (pup] 

+ztf[cos 6'df cos /3df sin (/)df + sin/3DF cos (/)df] (109) 

63 = Xtf sin 6'df sin /3df — Ytf [cos ^df sin /3df cos (pDp + cos ^df sin (^df] 

— ztf[cos 6'df sin /3df sin 0df — cos /3df cos (/)df] (1 10) 

or, equivalently: 

Xtf = ei cos 6'df — ^2 sin 6'df cos /3df + ^3 sin 6'df sin /3df (1 1 1) 

yTF = ei sin^DF cos (/)DF + e2[cos 0DF cos /3df cos ^DF — sin/3DF sin (/)df] 

—63 [cos 6'df sin /3df cos (^DF + cos /?df sin (/^df] (112) 
ztf = ei sin 0df sin (puF + ^2 [cos 6'df cos /3df sin 0df + sin /3df cos ^df] 

—63 [cos 6'df sin /3df sin 0df — cos /3df cos 0df] (113) 

Define the rotation matrix Rij = Xi ■ ej : 

R^J = (114) 

(COS &DF sin 6df cos (jiuF sin 6uf sin (j)uF \ 

-sin 0DF cos /3dp cos 6df cos/?dp cos (^df — sin /3df sin <^df cos 6df cos/3df sin ^oF + sin/^oF cos (^df | 
sin 9uF sin /?df —cos ^f sin /3df cos (/>df — cos /?df sin ^df —cos Oof sin /3df sin ^df + cos /3df cos ^df/ 

(115) 

and its inverse 

{R~'h= (116) 

(cos ^DF — sin ^F COS /9df sin ^df sin /3df \ 

sin ^DF cos <^DF COS 9f>f cos /3df COS (puF — sin /3df sin 4>f>f — cos 6'df sin /3df cos <^df — cos /3df sin 4>-df | 
sin^DF sin (/f>DF cos Sdf cos /3df sin c/iDF+sin /3df cos ^df — cos ^df sin /3df sin 0df + cos /3df cos 4>uf/ 

(117) 

The dielectric tensor ct is diagonal in the DF. If we calculate the dipole polarizabilility tensor in the DF 
it will also be diagonal, with elements 

( \ 

a°^= a^i (118) 
V a^i ) 

The polarizability tensor in the TF is given by 

(aTF),„ = i?.,(a°^),fe(i?-i)fc™ (119) 

Thus, we can describe a general anisotropic material if we provide, for each lattice site, the three 
diagonal elements of the dielectric tensor in the DF, and the three rotation angles 6'df, Pt)F, and 
<^DF- We first use the LDR prescription and the elements to obtain the three diagonal elements a^^ of 



the polarizability tensor in the DF. We then calculate the polarizability tensor a^^ using eq. (1119 1 
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28 Near-Field Calculations: Calculating E In or Near the Target 
28.1 Running DDSCAT 7.2 with NRFLD = 1 

For some scientific applications (e.g., surface-enhanced Raman scattering) one wishes to calculate the 
electric field E near the the target surface. It may also be of interest to calculate the electromagnetic field 
at positions within the target volume. DDSCAT 7.2 includes the capability for fast calculations of E in 
and near the target, using methods described by Flatau & Draine|p012) 



When DDSCAT is run with NRFLD=1, DDSCAT will be run twice. The first run, using a "minimal" 
computational volume just enclosing the physical target, creates stored files 

vxxxryyy'kzzz . poln 

for n=l and (if I0RTH=2) n=2. These files contain the stored solution for the polarization field Pj for 
lattice points within the minimal computational volume. 

DDSCAT 7.2 then takes the stored solution Pj and proceeds to calculate E at a lattice of points in 
an "extended" rectangular volume that can be specified to be larger than the original "minimal" compu- 
tational volume. The calculation is done rapidly using FFT methods, as described by |Flatau & Draine] 
( |2012| l. The result is then stored in output files 

Mxxxryyykzzz . En 
for n=l and (if I0RTH=2) n=2. 



28.2 The Binary Files wxxxryyykzzz . En 

For each of the points j — 1, N^yz in the extended volume where E was to be calculated, the binary 
file vxxxryyykzzz . E« contains: 

• The incident field Ejnc.j at each point j. 

• The field Egcaj at each point j produced by the polarization of the target (not including the dipole 
at j). (The total field at j is just E^ = Einc,j + EgcajO- 

• The polarization P j . Pj will be zero at all points outside the target (the ambient medium, taken to 
be vacuum). 

• The diagonal elements of the polarizability tensor aj at each point. An exact solution would 
satisfy Pj — otj (Eincj + ^sca.j)- Thus the stored P, a, Ejnc, and Esca allow the accuracy of the 
numerical solution to be verified. 

• The composition identifer /comp j at each lattice site. (The vacuum sites have /comp.j = 0). 

The vxxxryyykzzz . En files can be quite large. For example, the files wOOOrOOOkOOO.El and w000r000k000.E2 
created in the sample calculation in examples_exp/ELLIPSOID_NEARFIELD (see ^20.12.2 1 are 



each 86 MBbytes. Some of the stored data is easily recomputed (e.g., Ejnc and a.) or in principle re- 
dundant (Pj could in principle be obtained from ex, Ejnc, and Egca) but it is convenient to have them at 
hand. 



28.3 The program readnf 

A separate program, readnf . f 90, is provided to conveniently read the stored 

vxxxryyylizzz . En 

files. To create the readnf executable, position yourself in the /src directory and type 

make readnf 

which will compile READNF . f 90 and create an executable readnf. 
The program readnf is constructed to: 

1. Read a user-specified binary file (of the form wxxxryyykzzz . En) . The filename is specified in a 
parameter file readnf .par 

2. Read in control parameter I LINE (0 or 1) determining whether to calculate and write out E at 
points along a user-specified line. 
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3. Read in control parameter IVTR (0 or 1) determining whether to create files for subsequent visu- 
alization using VTK-based software. 

4. If ILINE=1 

• Read in seven parameters definining points along a line: XA,yA,ZA,XB,yB,ZB and Nab- 

• Obtain the complex E field at Nab uniformly- spaced points along the line connecting points 

{XA,yA,ZA) and {XB,yB,ZB)- 

• repeat for as many lines specifying xa, Va, za, XB,yB, zb and Nab as are present in the 
parameter file readnf . par. 

The program readnf writes output to an ascii output file readnf . out. The first 10 lines include 
information about the calculation, and column headings. Beginning with line 13, readnf . out gives 
physical coordinates x^p, y^p, ztf (in the "Target Frame"), and the real and imaginary parts of the x, 
y, and z components of the electric field E = Ei„c -|- E^ca at Nab points running from (a;^, yA,ZA) to 

{XB,yB,ZB)- 

We have carried out nearfield calculations for the problem of an Au sphere with radius a = 0.39789/xm 
(D = 0.7958/um) illuminated by a plane wave with A = 0.5/Ltm The Au target has refractive index m = 
0.96 + 1. Oil. The nearfield calculation of E is for a volume extending Q.5D beyond the spherical target in 
the +x,-x,+y,-y,+z,-z directions. This is the example problem in examples_exp/ELLIPSOID_NEARFIELD, 
where you can find the ddscat . par file. 

For the (x, 0. 0) track running through the center of the sphere the output readnf . out file wiU 
look like (for brevity, here we show only every 10th line of output along the track...): 

3.9789E-01 = a_eff {radius of equal -volume sphere) 
1.6408E-02 = d = dipole spacing 

59728 = N = number of dipoles in target or TUC 
96 96 96 = NX,NY,NZ = extent of computational volume 

5.0000E-01 = wavelength in vacuo 

1.00000 = refractive index of ambient medium 
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In the example given, the incident wave is propagating in the +x direction. The first layer of dipoles 
in the rectangular target is at x/d = —31.5, and the last layer is at x/d — —0.5. The "surface" of the 
target is at x/d = —32 (x = — 0.50/ini) and x/d = (x — 0). 



Figure 13 shows how \E\ varies along a line passing through the center of the Au sphere. As 



expected, the wave is strongly suppressed in the interior of the Au sphere. 
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Figure 13: Solid line: \E\'^ on track passing through center of Au sphere for the sample problem in exam- 
ples_exp/ELLIPSOID_NEARFIELD, calculated using target option ELLIPSOID. Note how weak the E field is inside the Au. 
Dashed line: \E\^ on track passing near the corner of the slab. Incident radiation is propagating in the +x direction. 



29 Displaying Target Shapes 

29.1 VTRCONVERT 

It is often desirable to be able to display the target shape. We are providing a program VTRCON- 
VERT.f90 which allows to convert DDSCAT shape format to VTK format. VTK is used world-wide in 
many advanced visualization applications such as: Para View, Visit, SDSlicer, or MayaVi2. 

Every time DDSCAT 7.2 is run, it will create a "target.out" file. To obtain a target.out file without 
running DDSCAT 7.2 the user can run calltarget to create a "target.out" file. For example 

calltarget < sphere40x40x40.shp 
where file "sphere40x40x40.shp" is 

ELLIPSOID 

40 40 40 

000 

will create an ASCII file "target.out" which is a list of the occupied lattice sites. The format of "target.out" 
is the same as the format of the "shape.dat" files read by DDSCAT if option FROM_FILE is used in 
ddscat.par (when you run the DDSCAT code). The DDSCAT format is very simple but is not compatible 
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with modern graphics programs such as "Para View" or "Mayavi2". Therefore we have created a program 
VTRCONVERT to convert between DDSCAT format and VTK format. 

The VTRCONVERT.f90 reads target shape data "target.out" and converts it to VTR format. The 
calUng sequence is 

VTRCONVERT target.out output 

where target.out is the name of the DDSCAT shape file created by "calltarget". The code writes 
two output files: "output_l.vtr" and "output.pvd". These files can be directly read in by "Para View" or 
"Mayavi2". 

If you are familar with PERL you can execute a script "shapes.pl". It will create several shape files. 

29.2 What is VTK? 

The Visualization Toolkit (VTK) is an open-source, freely available software system for 3D computer 
graphics and visualization. VTK is used world-wide in many advanced visualization applications such 
as: Para View, Visit, 3DSlicer, or MayaVi2. See for example 
http: // en . wikipedi a. org/ wiki/ParaView| 
'http://en.wikipedia.org/wiki/MayaVr 
Our converter code VTRCONVERT. f90 is written in FORTRAN90 and rehes on the pubhc domain 
software module vtr.f90 written by Jalel Chergui - see 

http://www.limsi.fr/Individu/chergui/pv/PVD.htm' 
vtrf90 defines 5 subroutines which enable writing ASCII data in XML VTK format. The code allows 
the user to write a 3D rectihnear mesh defined by x, y and z components. 

29.3 How to plot shapes once you have VTR/PVD files. 

Go to 'http://paraview.org/paraview/resources/software.html and download the Para View executable to 
your system (Para View is available for 'Windows, Linux, and MacOS) 
To plot a contour of a surface: 

1. In the toolbar, go to File/Open and select the file, e.g., "shapes/cylinder80x40.pvd " (note that in 
Para View one needs to open *.pvd file) 

2. Under Object Inspector/Properties click "apply" 
To add spheres indicating dipole positions 

1 . In the toolbar, go to Filters/Common/Contour 

2. under Object Inspector/Properties click "apply" 

3. In the toolbar, go to Filters/Common/Glyph 

4. Under Object Inspector/Properties one needs to change several settings 

5. Glyph Type - change to "sphere" 

6. Then 

(a) change "Radius" to 0.1 (make radius of a dipole smaller) as an initial choice. 

(b) change "Maximum Number of Points" to 5000 (or perhaps more depending on shape, but the 
points you have the longer it takes to plot) 

(c) toggle "Mask Points" to off 

(d) toggle "Random Mode" to off (plot dipoles in their positions) 

7. Click "Apply" in "object inspector" 

The initial choice of 0.1 for the radius in step 6 may not produce "touching" spheres - you can 
return to step 6 and iterate until you like the appearanace 

8. To rotate the figure, left-click on it and "drag" left-right and/or up-down to rotate around the vertical 
and/or horizontal axes. 

9. To add arrow showing target axes ai 



VISUALIZATION OF THE ELECTRIC HELD 



(a) In the toolbar, go to File/Open/ala2_l.pvd 

This file has just one point (at position (0,0,0)) and two vectors al(3) and a2(3) which are 
"target axes". 

(b) Under Object Inspector/Properties click "Apply". 

(c) In the Toolbar, go to Filters/Common/Glyph 

(d) Under Object Inspector/Properties 

i. "Vectors" will show "a 1" 

ii. change Glyph Type to "arrow" 

iii. change "Glyph Type" to "arrow" 

iv. change "Tip Radius" to 0.03 (improves appearance; suit yourself) 
V. change "Tip Length" to 0.1 (improves appearance...) 

vi. change "Shaft Radius" to 0.01 (improves appearance...) 

vii. cHck "edit" next to "Set Scale Factor" 

viii. change "Scale Factor" to some value like "40" 
ix. toggle "Mask Points" to off 

X. toggle "Random Mode" to off 
xi. click "Apply" 

10. To add arrow showing target axis a2 

(a) In the toolbar, go to File/Ope n/ala2_l.pvd 

(b) Under Object Inspector/Properties click "Apply". 

(c) In the Toolbar, go to Filters/Common/Glyph 

(d) Under Object Inspector/Properties 

i. "Vectors" will show "al" - change this to "a2" 

ii. change "Glyph Type" to "arrow" 

iii. change "Tip Radius" to 0.03 (improves appearance...) 

iv. change "Tip Length" to 0.1 (improves appearance...) 

V. change "Shaft Radius" to 0.01 (improves appearance...) 

vi. click "edit" next to "Set Scale Factor" 

vii. change "Scale Factor" to some value like "40" 

viii. toggle "Mask Points" to off 
ix. toggle "Random Mode" to off 
X. click "Apply" 

11. To save the image: in the toolbar, go to File/Save Screenshot 

(a) select desired resolution (default may be acceptable) and click "Ok". 

(b) enter the desired filename 

(c) select the file type (options are .jpg, .png, .pdf, .tif, .ppm, .bmp) 

(d) click "OK" 

If you wish to add additional features to images, please consult the Para View documentation. 



30 Visualization of the Electric Field 

The following instructions assume the user has the MayaVi2 graphics package installed. If you don't 
have it yet, go to 

http://code.enthoug ht.com/projects/mayavi/mayavi/installation.htmll 
for installation options. 

1. First mn example ELLIPSOID_NEARFIELD. This will produce output_l.vtr file. 
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Figure 14: Visualization witli Para View of tlie dipole realization of a sphere, produced following above instructions, using files 
in directory examples_exp/ELLIPSOID. 

2. Start Mayavi2 

3. File/load data/open file output 

4. Point to examples_exp/ELLlPSlD_NEARFlELDA^TRoutput_l.vtr 

5. colors and legends/add module/outline 

6. outline/right click/add module/contour grid plane 

7. In Mayavi2 object editor you can position slider between (0,95). Choose value in the middle. Set 
"filled contours" to on. Increase number of contours to 64. 

8. One can add more "contour grid planes" to illustrate the results in different cross sections. 

31 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the 
ddscat . par file is intended to be simple and suggestive so that, after reading the above notes once, 
the user may not have to refer to them again. 

Known bugs in DDSCAT will be posted at the DDSCAT web site, 

http : / / www . astro . princeton . edu / --'draine/ DDSCAT . html 
and the latest version of DDSCAT will always be available at this site. The source code is also available 
at 

|http : / /code . google . com/p/ddscat /| 
There is also a wiki: 

|http : / /ddscat . wikidot . com/start| 

Users are encouraged to provide B. T. Draine (draine@astro.princeton.edu) with their 
email address; email notification of bug fixes, and any new releases of DDSCAT, will be made known 
to those who do! 

P. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" with URL 
|http://at ol . ucsd . edu/~^pf latau The SCATTERLIB Internet site is a library of light scat- 
tering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information 
related to scattering on spherical and non-spherical particles is collected: an extensive list of references 
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Figure 15: |E|/|Eo| on two planes, one passing through the center and one passing near the a — 0.398/im Au sphere calculated 
in examples_exp/ELLIPSOID_NEARFIELD. The incident wave, with A = 0.5fim, is propagating with ko || xtf- and Eo || ytf- 
The axes are labelled in /im. This is the same problem as the results shown in FigurefTS] This figure was generated by MayaVi2. 



to light scattering methods, refractive index, etc. This URL page contains section on the discrete dipole 
approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. 

Users of DDSCAT should cite appropriate papers describing DDA theory and its implementation. 
The following papers may be relevant (pdfs are included in the /doc directory) 



Draine ( 1988| l: basic theory of the DDA, including radiative reaction, and application to anisotropic 



materials; 



Goodman et al. ( 1990 1: introduction of FFT methods to greatly accelerate DDA calcualations 

[Draine & Flatau| ( |2004[ ): review of the DDA, including demonstrations of accuracy and conver- 
gence; 



Draine & Flatau ( 2008| l: extension of the DDA to periodic structures, and generalization of the 



Mueller scattering matrix to describe scattering by 1-d and 2-d periodic structures. 
Flatau & Draine ( 2012[ l: implementation of efficient near-field calculations using FFTs. 



• This UserGuide, for which we suggest the following citation: 

Draine, B.T., & Flatau, P.J. 2012, "User Guide for the Discrete Dipole Approximation Code 



DDSCAT (Version 7.2)' Mhttp://arxiv.org/abs/ 1202.3424 



Finally, the authors have one special request: We would appreciate receiving copies (pdf is fine) of 
papers which make use of DDSCAT! 
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A Understanding and Modifying dd scat. par 

In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to 
modify the ddscat . par file. Here we list the sample ddscat . par file for the example problem in 
examples_exp/RCTGLPRSM, followed by a discussion of how to modify this file as needed. Note 
that all numerical input data in DDSCAT is read with free-format READ ( IDEV, * ) . . . statements. 
Therefore you do not need to worry about the precise format in which integer or floating point numbers 
are entered on a line. The crucial thing is that fines in ddscat . par containing numerical data have the 
correct number of data entries, with any informational conoments appearing after the numerical data on 
a given fine. 

' ========= Parameter file for v7 . 2 ===================' 

Preliminaries 

'NOTORQ' = CMT0RQ*6 (NOTORQ, DOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, GPBICG, PETRKP, QMRCCG) — solution method 
'GPFAFT' = CMDFFT*6 (GPFAFT, FFTMKL) — FFT method 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) — prescription for polarizabilities 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) — specify binary output 
' Initial Memory Allocation 

100 100 100 = dimensioning allowance for target generation 
' **** Target Geometry and Composition 
'RCTGLPRSM' = CSHAPE*9 shape directive 
16 32 32 = shape parameters 1-3 

1 = NCOMP = number of dielectric materials 
' . . /diel/Au_evap' = file with refractive index 1 

Additional Nearfield calculation? 

= NRFLD (=0 to skip nearfield calc, =1 to calculate nearfield E) 

0.0 0.0 0.0 0.0 0.0 0.0 (fract. extens. of calc. vol. in -x, +x, -y, +y, -z, +z) 
Error Tolerance 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC| E>-ACA | X> ) / (NORM OF AC|E>) 
' -k-k-k-k maximum number of iterations allowed 
300 = MXITER 

' **** Interaction cutoff parameter for PBC calculations 

l.OOe-2 = GAMMA (le-2 is normal, 3e~3 for greater accuracy) 

' ±kkk Angular resolution for calculation of <cos>, etc. 

0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] "2 ) 

Vacuum wavelengths (micron) 
0.5000 0.5000 1 'LIN' = wavelengths (first, last, how many, how=LIN, INV, LOG) 
' kkk* Refractive index of ambient medium' 
1.000 = NAMBIENT 

'**** Effective Radii (micron) **** ' 

0.246186 0.246186 1 'LIN' = aeff (first, last, how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol. state eOl; =2 to also do orth. pol . state) 
' *kkk Specify which output files to write ****' 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

'-k-k-k-k Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 0. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'k-k-k-k Select Elements of S_ij Matrix to Print kkkk' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 ^ indices ij of elements to print 
'-kkkk Specify Scattered Directions ****' 

'LFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

2 = NPLANES = number of scattering planes 

0. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 1 
90. 0. 180. 5 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 2 
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Lines Comments 
1-2 comment lines 

3 NOTORQ if torque calculation is not required; 

DOTORQ if torque calculation is required. 



4 PBCGS2 is recommended; other options are PBCGST, GPBICG, PETRKP, and QMRCCG (see { 12 1. 

5 GPFAFT is suppUed as default, but FF TMK L is recommended if DDSCAT has been compiled with 

)n lot 



the Intel® Math Kernel Library (see §S 6.5 



6 GKDLDR is recommended as the prescription for polarizabilities (see 5 14 1 

7 NOTE IN for no unformatted binary output. 

ORIBIN for unformatted binary dump of ori entati onal averages only; 
ALLBIN for full unformatted binary dump (5 1 1.2 >; 

8 comment line 

9 initial memory allocation NX,NY,NZ. These must be large enough to accomodate the target that 
will be generated. 

10 comment line 

1 1 specify choice of target shape (see 5 20 for description of options RCTGLPRSM, ELLIPSOID, 
TETRAHDRN, ...) 

12 shape parameters SHPARl, SHPAR2, SHPAR3, ... (see d20 . 

13 number of different dielectric constant tables (see S 15 i. 

14 name(s) of dielectric constant table(s) (one per line). 

15 comment line 

16 NRFLD = (or 1) to skip (or do) nearfield calculation 

17 6 non-negative numbers ri, ra specifying fractional extensi on o f comptutational volume 
(in — xtf, +xtf, — Ytf, +yTF, — ztf, +ztf direction (see i 9.4 1 

18 commment line 

20 TOL = error tolerance h: maximum allowed value of \A'' E — A'' AP\/\A'' E\ [see eq.jTsj]. 

21 comment line 

22 MXITER= maximum number of conjugate-gradient iterations allowed 

23 comment line 

24 GAMMA= interaction cutoff parameter 7 (see Draine & Flatau 2009) 
the value of 7 does not affect calculations for isolated targets 

25 comment line 

26 ETASCA - parameter ri controlling angular averages (i 24 1. 

27 comment line 

28 A - vacuum wavlenghts: first, last, how many, how chosen. 

29 comment line 

30 NAMBIENT= (real) refractive index of ambient medium. 

3 1 comment line 

32 flcff - first, last, how many, how chosen. 

33 comment line 

34 specify x,y,z components of (complex) incident polarization eoi (^23[> 

35 lORTH = 2 to do both polarization states (normal); 
lORTH = 1 to do only one incident polarization. 

36 comment line 

37 IWRKSC = to suppress writing of ".sea" files; 
IWRKSC = 1 to enable writing of ".sea" files. 

38 comment line 

39 /3 (see 5 Tsl - first, last, how many . 

40 G (see i - first, last, how many. 

41 $ (see 5 Tsl - first, last, how many. 

42 comment line 

43 IWAVO IRADO lORIO ~ Starting values of integers IWAV IRAD I OR I (normally 0). 

44 comment line 

45 Ns = number of scattering matrix elements (must be < 9) 

46 indices ij of Ns elements of the scattering matrix Sij 

47 comment line 

48 specify whether scattered directions are to be specified by CMDFRM=' LFRAME' or ' TFRAME' . 

49 NPLANES= number of scattering planes to follow 

50 (pa for first scattering plane, Os.min, Os,max, how many 8s values; 
51,... (j)s for 2nd,... scattering plane, ... 
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B vrxxxryyy . avg Files 



The file wOOOrOOOori. avg contains the results for the first wavelength (wOOO) and first target radius 
(rOOO) averaged over orientations (.avg). The wOOOrOOOori . avg file generated by the sample 
calculation in examples_exp/RCTGLPRSM should look Uke the following: 

DDSCAT DDSCAT 7.2.0 [12.02.11] 

TARGET Rectangular prism; NX,NY,NZ= 16 32 32 

PBCGS2 method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

^ NATO ^ number of dipoles 

= d/aeff for this target [d=dipole spacing] 

= d (physical units) 

0.246185 ^ effective radius (physical units) 
0.500000 ^ wavelength (in vacuo, physical units) 
3.093665 = 2*pi*aeff /lambda 

1.000000 = refractive index of ambient medium 

1.8628), eps.= ( -2.5374 , 3.5975) |m|kd= 0.4120 for subs. 1 
OOOE-05 = error tolerance for CCG method 
0.00000 0.00000) = target axis Al in Target Frame 
1.00000 0.00000) = target axis A2 in Target Frame 
962 = (theta,phi) values used in comp. of Qsca,g 
0.00000 0.00000) = k vector (latt. units) in Lab Frame 
0.00000) ( 1.00000, 0.00000) ( 0.00000, 0.00000)=inc.pol.vec. 1 in LF 
0.00000, 0.00000) (-0.00000, 0.00000) ( 1. 00000, -O.00000)=inc.pol.vec. 2 in LF 
0.000 0.000 = beta_min, beta_max ; NBETA = 1 
0.000 = theta_min, theta_max; NTHETA= 1 
0.000 = phi_min, phi_max ; NPHI = 1 



16384 
.06346820 
.015625 
AEFF= 
WAVE= 
K*AEFF= 
NAMBIENT= 
n= ( 0.9656 

TOL= 1 
( 1.00000 
( 0.00000 

NAVG= 
( 0.19635 
( 0.00000 
( 

.000 
.000 



0.5000 = ETASCA = param. controlling # of scatt. dirs used to calculate <cos> etc. 
Results averaged over 1 target orientations 

and 2 incident polarizations 
Qext Qabs Qsca g(l)=<cos> <cos'^2> Qbk Qpha 

J0=1: 3.6135E+00 1.4308E+00 2.1827E+00 4.1463E-01 7.5169E-01 6.1831E-01 -1.2614E-01 
J0=2: 3.6135E+00 1.4308E+00 2.1827E+00 4.1463E-01 7.5169E-01 6.1831E-01 -1.2614E-01 
mean: 3.6135E+00 1.4308E+00 2.1827E+00 4.1463E-01 7.5169E-01 6.1831E-01 -1.2614E-01 
Qpol= -1.6212E-05 dQpha= 1.2815E-06 

Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter Nsca 
J0=1: 9.0502E-01 -1.0545E-08 1.6702E-06 18 300 962 

J0=2: 9.0503E-01 7.8238E-07 9.9115E-08 18 300 962 

mean: 9.0503E-01 3.8592E-07 8.8467E-07 

Mueller matrix elements for selected scattering directions in Lab Frame 
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-06 


80, 


.00 


90, 


.00 


0. 


.83815 


1 , 


.4711E+00 


1 , 


.2330E+00 


1. 


.233E+00 


1 . 


.471E+00 


8. 


. 662E- 


-06 


5. 


. 942E- 


-06 


85, 


.00 


90, 


.00 


0. 


.87905 


1 , 


.3178E+00 


1 , 


.1584E+00 


1. 


.158E+00 


1. 


.318E+00 


8. 


.342E- 


-06 


6. 


.288E- 


-06 


90, 


.00 


90, 


.00 


0. 


.84389 


1 , 


.1878E+00 


1 , 


. 0024E+00 


1. 


. 002E+00 


1. 


.188E+00 


7. 


. 830E- 


-06 


6. 


.248E- 


-06 


95, 


.00 


90, 


.00 


0. 


.78582 


1 , 


.0596E+00 


8, 


.3269E-01 


8. 


. 327E-01 


1 . 


.060E+00 


7 . 


.320E- 


-06 


6. 


. 148E- 


-06 


100, 


.00 


90 , 


.00 


0. 


. 75673 


9, 


.2364E-01 


6, 


. 9895E-01 


6. 


. 989E-01 


9. 


. 236E-01 


6. 


. 717E- 


-06 


5 . 


. 727E- 


-06 


105 , 


.00 


90 , 


.00 


0. 


.78540 


7 , 


. 9129E-01 


6, 


.2148E-01 


6. 


. 215E-01 


7 . 


. 913E-01 


6. 


. 345E- 


-06 


5 . 


.222E- 


-06 


110, 


.00 


90 , 


. 00 


. 


.85741 


6 , 


. 9404E-01 


5 , 


. 9507E-01 


5. 


. 951E-01 


6. 


. 940E-01 


5 . 


. 973E- 


-06 


4 . 


. 621E- 


-06 


115 , 


.00 


90 , 


.00 


0. 


.89355 


6, 


. 7419E-01 


6, 


. 0242E-01 


6. 


. 024E-01 


6. 


. 742E-01 


5 . 


. 681E- 


-06 


3. 


. 962E- 


-06 


120, 


.00 


90, 


.00 


0. 


.81395 


7 , 


. 7423E-01 


6, 


. 3018E-01 


6. 


. 302E-01 


7 . 


. 742E-01 


5 . 


.308E- 


-06 


3. 


.534E- 


-06 


125, 


.00 


90, 


.00 


0. 


. 65812 


1 , 


. 0320E+00 


6, 


. 7920E-01 


6. 


. 792E-01 


1 . 


. 032E+00 


4 . 


. 740E- 


-06 


3. 


. 397E- 


-06 


130 , 


.00 


90, 


.00 


0. 


. 51335 


1 , 


. 4856E+00 


7 , 


. 6265E-01 


7 . 


. 627E-01 


1 . 


. 486E + 00 


3. 


. 845E- 


-06 


4 . 


. 038E- 


-06 


135 , 


.00 


90, 


.00 


0. 


.40829 


2 , 


. 1850E+00 


8, 


. 9213E-01 


8. 


. 921E-01 


2 . 


. 185E+00 


2 . 


.518E- 


-06 


5 . 


.262E- 


-06 


140 , 


.00 


90, 


.00 


0. 


.33090 


3, 


. 2005E+00 


1 , 


. 0590E+00 


1 . 


. 059E + 00 


3. 


. 200E+00 


6. 


. 786E- 


-07 


7 . 


.379E- 


-06 


145 , 


.00 


90, 


.00 


0. 


.26513 


4 , 


. 6140E+00 


1 , 


. 2233E+00 


1 . 


. 223E + 00 


4 . 


. 614E + 00 


-1 . 


. 627E- 


-06 


1 . 


. 036E- 


-05 


150 , 


.00 


90, 


.00 


0. 


.20351 


6, 


. 4859E+00 


1 , 


. 3199E+00 


1 . 


. 320E+00 


6. 


. 486E + 00 


-4 . 


. 127E- 


-06 


1 . 


. 380E- 


-05 


155, 


.00 


90, 


.00 


0. 


. 14592 


8, 


. 8035E+00 


1 , 


.2846E+00 


1 . 


. 285E + 00 


8. 


. 804E + 00 


-6. 


. 745E- 


-06 


1 . 


. 739E- 


-05 


160, 


.00 


90, 


.00 


0. 


.09522 


1 , 


.1432E+01 


1 , 


.0885E+00 


1. 


.089E+00 


1 . 


.143E+01 


-9. 


. 197E- 


-06 


2 . 


. 079E- 


-05 


165, 


.00 


90, 


.00 


0. 


.05409 


1 , 


.4098E+01 


7 , 


. 6255E-01 


7 . 


. 625E-01 


1. 


.410E+01 


-1. 


. 131E- 


-05 


2 . 


.337E- 


-05 


170, 


.00 


90, 


.00 


0. 


. 02413 


1 , 


. 6425E+01 


3, 


.9640E-01 


3. 


. 964E-01 


1. 


. 642E + 01 


-1. 


.240E- 


-05 


2 . 


.484E- 


-05 


175, 


.00 


90, 


.00 


0. 


.00604 


1 , 


.8022E+01 


1 , 


.0878E-01 


1 . 


.088E-01 


1 . 


.802E+01 


-1. 


.315E- 


-05 


2 . 


.488E- 


-05 


180, 


.00 


90, 


.00 


0. 


.00000 


1 , 


.8591E+01 


-6, 


.6757E-06 


-6. 


.676E-06 


1 . 


.859E+01 


-1. 


.306E- 


-05 


2 . 


.314E- 


-05 
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C wxxxryyykzzZ'Sca Files 

The wOOOrOOOkOOO . sea file contains the results for the first wavelength (wOOO), first target radius 
(rOOO), and first orientation (kOOO). The wOOOrOOOkOOO . sea file created by the sample calculation 
in examples_exp/RCTGLPRSM should look like the following: 

DDSCAT DDSCAT 7.2.0 [12.02.11] 

TARGET Rectangular prism; NX,NY,NZ= 16 32 32 

PBCGS2 method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

16384 = NATO = number of dipoles 

0.06346820 = d/aeff for this target [d=dipole spacing] 
0.015625 = d (physical units) 

physical extent of target volume in Target Frame 

-0.250000 0.000000 = xmin,xmax (physical units) 

-0.250000 0.250000 = ymin,ymax (physical units) 

-0.250000 0.250000 = zmin, zmax (physical units) 

0.246186 = effective radius (physical units) 
0.500000 = wavelength (in vacuo, physical units) 
3.093665 = 2*pi*aeff /lambda 

1.000000 = refractive index of ambient medium 

1.8628), eps.= ( -2.5374 , 3.5975) |m|kd= 0.4120 for subs. 1 
OOOE-05 = error tolerance for CCG method 
0.00000 0.00000) = target axis Al in Target Frame 



AEFF= 

WAVE= 
K*AEFF= 
NAMBIENT= 
n= ( 0.9656 

TOL= 1 
( 1.00000 
( 0.00000 

NAVG= 



1.00000 0.00000) = target axis A2 in Target Frame 
962 = (theta,phi) values used in comp. of Qsca,g 
( 0.19635 0.00000 0.00000) = k vector (latt . units) in TF 
( 0. 00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 
( 0. 00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, O.00000)=inc.pol.vec. 
( 1.00000 0.00000 0.00000) = target axis Al in Lab Frame 
(-0.00000 1.00000 0.00000) = target axis A2 in Lab Frame 
( 0.19635 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0. 00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 
( 0. 00000, 0.00000) (-0.00000, 0.00000) ( 1. 00000, -O.00000)=inc.pol.vec. 
0.000 = rotation of target around Al 
0.000 = angle between Al and k 
0.000 = rotation of Al around k 



1 in TF 

2 in TF 



1 in LF 

2 in LF 



BETA = 
THETA= 
PHI = 
0.5000 



J0=1 
J0=2 
mean 
Qpol^ 

J0=1 
J0=2 
mean 



= ETASCA = 

Qext 
3. 6135E+00 
3. 6135E+00 
3. 6135E+00 
1 . 6212E-05 
Qsca*g ( 1 ) 
9.0502E-01 
9.0503E-01 
9.0503E-01 



param. controlling # of scatt . dirs used to calculate <cos> etc. 



Qabs 

L.4308E+00 2 

L.4308E+00 2 

L.4308E+00 2 

Qsca*g (2) 
-1.0545E-08 
7.8238E-07 
3.8592E-07 



Qsca 

. 1827E+00 
.1827E+00 4 
.1827E+00 4 

Qsca*g ( 3 ) 
1 . 6702E-06 
9.9115E-08 
8.8467E-07 



g ( 1 ) =<cos> 
4 . 1463E-01 



<cos^2> 
7 . 5169E-01 



Qbk 
6.1831E-01 



1463E-01 
1463E-01 



7.5169E-01 6.1831E-01 
7.5169E-01 6.1831E-01 
dQpha= 



Qpha 
-1 . 2614E 
-1.2614E 
-1.2614E 

1.2815E 



-01 
-01 
01 
06 



iter 
18 
18 



mxiter 
300 
300 



Nsca 
962 
962 



Mueller matrix elements for selected scattering directions in Lab Frame 



theta 


phi 




Pol . 


S_ll 


S_12 


S_21 


S_22 




3_31 






3_41 







00 





.00 





00000 


7 


.5116E+01 


3 


4332E-04 


3 


433E-04 


7 


.512E+01 


1 


. 646E- 


05 


-2 


. 015E- 


05 


5 


00 





.00 





00767 


7 


.2167E+01 


5 


5379E-01 


5 


538E-01 


7 


.217E+01 


2 


. 202E- 


05 


-2 


. 682E- 


05 


10 


00 





.00 





03116 


6 


.3969E+01 


1 


9930E+00 


1 


993E+00 


6 


.397E+01 


1 


. 927E- 


05 


-3 


. 127E- 


05 


15 


00 





.00 





07180 


5 


.2244E+01 


3 


7511E+00 


3 


751E+00 


5 


.224E+01 


1 


.716E- 


05 


-3 


. 079E- 


05 


20 


00 





00 





13159 


3 


9245E+01 


5 


1644E+00 


5 


164E+00 


3 


. 924E+01 


1 


. 876E- 


05 


-2 


. 890E- 


05 


25 


00 





00 





21216 


2 


7088E+01 


5 


7468E+00 


5 


747E+00 


2 


709E+01 


1 


. 541E- 


05 


-2 


. 897E- 


05 


30 


00 





00 





31112 


1 


7234E+01 


5 


3617E+00 


5 


362E+00 


1 


723E+01 


9 


. 466E- 


06 


-2 


. 667E~ 


05 


35 


00 





00 





41156 


1 


0271E+01 


4 


2270E+00 


4 


227E+00 


1 


027E+01 


6 


.514E- 


06 


-2 


.221E- 


05 


40 


00 





00 





46174 


6 


0116E+00 


2 


7758E+00 


2 


776E+00 


6 


. 012E + 00 


4 


.550E- 


06 


-1 


. 549E- 


05 


45 


00 





00 





38394 


3 


7970E+00 


1 


4578E+00 


1 


458E+00 


3 


.797E+00 


2 


. 155E- 


06 


-1 


. 198E- 


05 


50 


00 





00 





20310 


2 


8383E+00 


5 


7646E-01 


5 


765E-01 


2 


. 838E+00 


3 


. 539E- 


06 


-8 


. 333E- 


06 


55 


00 





00 





08917 


2 


4759E+00 


2 


2078E-01 


2 


208E-01 


2 


.476E+00 


4 


.546E- 


06 


-8 


. 382E- 


06 


60 


00 





00 





12808 


2 


2936E+00 


2 


9378E-01 


2 


938E-01 


2 


294E+00 


5 


. 510E- 


06 


-6 


. 543E- 


06 


65 


00 





00 





28428 


2 


1074E+00 


5 


9909E-01 


5 


991E-01 


2 


107E+00 


7 


. 906E- 


06 


-7 


. 080E- 


06 


70 


00 





00 





49621 


1 


8862E+00 


9 


3596E-01 


9 


360E-01 


1 


.886E+00 


1 


. 034E- 


05 


-7 


.218E- 


06 


75 


00 





.00 





70066 


1 


. 6630E+00 


1 


1652E+00 


1 


165E+00 


1 


.663E+00 


1 


.263E- 


05 


-7 


.485E- 


-06 


80 


00 





.00 





83815 


1 


.4710E+00 


1 


2329E+00 


1 


233E+00 


1 


.471E+00 


1 


.378E- 


05 


-7 


.399E- 


06 


85 


00 





.00 





87905 


1 


.3178E+00 


1 


1584E+00 


1 


158E+00 


1 


.318E+00 


1 


.430E- 


05 


-7 


. 156E- 


06 
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89 



90, 


.00 


0, 


.00 


0. 


.84389 


1 , 


. 1878E+00 


1 . 


. 0024E + 00 


1 . 


. 002E+00 


1 . 


. 188E + 00 


1 . 


. 409E- 


-05 


-6. 


. 670E- 


-06 


95, 


.00 


, 


. 00 


0. 


.78583 


1 , 


. 0596E+00 


8. 


.3269E-01 


8 . 


. 327E-01 


1 . 


. 060E+00 


1 . 


.292E- 


-05 


-6. 


.226E- 


-06 


100, 


.00 





. 00 


. 


. 75672 


9, 


.2363E-01 


6. 


. 9893E-01 


6. 


. 989E-01 


9. 


. 236E-01 


1 . 


. 119E- 


-05 


-5 . 


. 588E- 


-06 


105. 


.00 


0, 


.00 


0. 


.78541 


7 , 


. 9129E-01 


6. 


.2148E-01 


6. 


.215E-01 


7 . 


. 913E-01 


1. 


. 052E- 


-05 


-4 . 


.689E- 


-06 


110. 


.00 


0, 


.00 


0. 


. 85741 


6, 


. 9404E-01 


5. 


. 9508E-01 


5. 


.951E-01 


6. 


. 940E-01 


9. 


. IIIE- 


-06 


-3. 


. 677E- 


-06 


115. 


.00 


0, 


.00 


0. 


.89355 


6, 


.7418E-01 


6. 


.0241E-01 


6. 


.024E-01 


6. 


.742E-01 


8. 


. 132E- 


-06 


-2 . 


.382E- 


-06 


120. 


.00 


0, 


.00 


0. 


.81394 


7 , 


.7423E-01 


6. 


.3018E-01 


6. 


.302E-01 


7 . 


.742E-01 


7. 


. 871E- 


-06 


-1 . 


.250E- 


-06 


125. 


.00 


0, 


.00 


0. 


.65812 


1 , 


.0320E+00 


6. 


.7920E-01 


6. 


.792E-01 


1 . 


.032E+00 


8. 


.355E- 


-06 


1 . 


. 059E- 


-06 


130. 


.00 


0, 


.00 


0. 


.51335 


1 , 


.4856E+00 


7 . 


. 6265E-01 


7 . 


. 626E-01 


1 . 


. 486E + 00 


9. 


.331E- 


-06 


3. 


.307E- 


-06 


135. 


.00 


0, 


.00 


0. 


.40829 


2 , 


.1850E+00 


8. 


. 9212E-01 


8. 


. 921E-01 


2 . 


.185E+00 


1. 


. 104E- 


-05 


7 . 


.216E- 


-06 


140. 


.00 


0, 


.00 


0. 


.33090 


3, 


.2005E+00 


1 . 


.0590E+00 


1. 


.059E+00 


3. 


.200E+00 


1. 


.369E- 


-05 


1 . 


. 023E- 


-05 


145. 


.00 


0, 


.00 


0. 


.26513 


4 , 


. 6140E+00 


1 . 


.2233E+00 


1. 


.223E+00 


4 . 


. 614E+00 


1. 


. 619E- 


-05 


1 . 


.302E- 


-05 


150. 


.00 


0, 


.00 


0. 


.20350 


6, 


.4859E+00 


1 . 


.3199E+00 


1 . 


.320E+00 


6. 


. 486E + 00 


1 . 


. 870E- 


-05 


1 . 


. 934E- 


-05 


155 , 


.00 


, 


.00 


0. 


.14592 


8 , 


. 8035E+00 


1 . 


. 2846E + 00 


1 . 


. 285E+00 


8 . 


. 804E + 00 


1 . 


. 867E- 


-05 


2 . 


.269E- 


-05 


160 , 


.00 


0, 


.00 


0. 


.09522 


1 , 


. 1432E+01 


1 . 


. 0885E+00 


1 . 


. 089E + 00 


1 . 


. 143E+01 


1 . 


. 872E- 


-05 


2 . 


. 696E- 


-05 


165 , 


.00 


, 


.00 


0. 


. 05409 


1 , 


. 4098E+01 


7 . 


. 6256E-01 


7 . 


. 626E-01 


1 . 


. 410E+01 


1 . 


. 602E- 


-05 


2 . 


. 858E- 


-05 


170 , 


.00 


, 


.00 


0. 


. 02413 


1 , 


. 6425E+01 


3. 


. 9637E-01 


3. 


. 964E-01 


1 . 


. 642E + 01 


1 . 


.593E- 


-05 


2 . 


. 845E- 


-05 


175 , 


.00 


0, 


.00 


0. 


. 00604 


1 , 


. 8022E+01 


1 . 


. 0884E-01 


1 . 


. OBBE-Ol 


1 . 


. 802E+01 


1 . 


. 408E- 


-05 


2 . 


. 652E- 


-05 


180 , 


.00 


0, 


.00 


0. 


.00000 


1 , 


. 8591E+01 


5 . 


. 7221E-06 


5. 


. 722E-06 


1 . 


. 859E+01 


1 . 


.300E- 


-05 


2 . 


. 317E- 


-05 


, 


.00 


90, 


.00 


0. 


.00000 


7 , 


. 5116E+01 


-3 . 


. 4332E-04 


-3. 


. 433E-04 


7 . 


. 512E+01 


-1 . 


. 684E- 


-05 


-2 . 


.006E- 


-05 


5 , 


.00 


90, 


.00 


0. 


.00768 


7 , 


.2167E+01 


5 . 


. 5393E-01 


5. 


. 539E-01 


7 . 


. 217E+01 


-1 . 


. 144E- 


-05 


-2 . 


. 152E- 


-05 


10, 


.00 


90, 


.00 


0. 


.03116 


6, 


. 3969E+01 


1 . 


. 9931E+00 


1 . 


. 993E + 00 


6 . 


. 397E + 01 


-6. 


. 057E- 


-06 


-2 . 


. 180E- 


-05 


15 , 


.00 


90, 


.00 


0. 


. 07180 


5, 


.2244E+01 


3. 


. 7512E+00 


3. 


. 751E+00 


5. 


.224E+01 


-1. 


.538E- 


-06 


-2 . 


. 022E- 


-05 


20. 


.00 


90, 


.00 


0. 


.13159 


3, 


. 9245E+01 


5. 


.1643E+00 


5. 


. 164E + 00 


3. 


. 924E + 01 


4 . 


.089E- 


-06 


-1 . 


. 916E- 


-05 


25. 


.00 


90, 


.00 


0. 


.21216 


2 , 


.7088E+01 


5. 


.7470E+00 


5. 


.747E+00 


2 . 


.709E+01 


6. 


. 62 6E- 


-06 


-1 . 


.703E- 


-05 


30. 


.00 


90, 


.00 


0. 


.31112 


1 , 


.7234E+01 


5. 


.3617E+00 


5. 


.362E+00 


1. 


.723E+01 


8. 


.549E- 


-06 


-1 . 


.415E- 


-05 


35. 


.00 


90, 


.00 


0. 


.41156 


1 , 


.0271E+01 


4 . 


.2270E+00 


4 . 


.227E+00 


1 . 


.027E+01 


9. 


.798E- 


-06 


-1 . 


. 066E- 


-05 


40. 


.00 


90, 


.00 


0. 


.46175 


6, 


.0116E+00 


2 . 


.7759E+00 


2 . 


.776E+00 


6. 


.012E+00 


9. 


. 635E- 


-06 


-7 . 


.258E- 


-06 


45. 


.00 


90, 


.00 


0. 


.38394 


3, 


.7970E+00 


1 . 


.4578E+00 


1. 


.458E+00 


3. 


.797E+00 


9. 


.749E- 


-06 


-4 . 


.346E- 


-06 


50. 


.00 


90, 


.00 


0. 


.20311 


2 , 


.8383E+00 


5. 


.7649E-01 


5. 


.765E-01 
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VlXXXRYYYKZZZ .POLN FILES 



D vjxxxryyykzzz .poln Files 

Binary files wxxxryyykzzz .poln are written to disk only when nearfield calculations are done (pa- 
rameter NRFLD=1). The wOOOrOOOkOOO .poll file contains the polarization solution for the first 
wavelength (wOOO), first target radius (rOOO), first orientation (kOOO), and first incident polarization 
(poll). In order to limit the size of this file, it has been written as an unformatted or "binary" file. 
This preserves full machine precision for the data, is quite compact, and can be read efficiently, but 
unfortunately the file is not fully portable because different computer architectures (e.g., Linux vs. MS 
Windows) have adopted different standards for storage of "unformatted" data. However, anticipating that 
many users will be computing within a single architecture, the distribution version of DDSCAT uses this 
format. 

Additional warning: even on a single architecture, users should be alert to the possibility that different 
compilers may follow different conventions for reading/writing unformatted files. 
The file contains the following information: 

• The location of each dipole in the target frame. 

• {k^, ky, kz)d, where k is the incident k vector 

• {Eox, Eoy, Eqz), the complex polarization vector of the incident wave. 

• a^^d'^, the inverse of the symmmetric complex polarizability tensor for each of the dipoles in the 
target. 

• {Px ,Py,Pz), the complex polarization vector for each of the dipoles. 

The interested user should consult the routine writepol . f to see how this information has been orga- 
nized in the unformatted file. 



E wxxxryyykzzz.En Files 

Binary files wxxxryyy]zzzz . En are written to disk only when nearfield calculations are done (parameter 
NRFLD=1). The wOOOrOOOkOOO .El file contains information describing the problem and the solution 
at "grid points" throughout the extended "computational volume" specified for the nearfield calculation 



(see 5 28 1. These binary files are very large, but have been written in a way to simplify subsequent use 



for visualization using, e.g., the program READNF . f 90 (see 5 30i. The interested user can examine 
subroutine nearfield . f 90 that writes the file, or program READNF . 190 that reads the file, to see 
how the information is organized. 

A user who finds the file size to be a serious problem may wish to modify the code in nearfield. 190 
to suppress writing of some of the arrays (e.g., the diagonal elements of the A matrix), retaining only 
the data of specific interest (e.g., the E field). Of course, and modifications to WRITE statements in 
nearfield .190 will require corresponding changes to READ statements in READNF . f 90. 
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Qabs, 11 
Qoxt, 1 1 
Qsca,^ 1 1 
<3pha. 11 

Sj - scattering amplittide matrix, 68 

Sij - 4x4 Mueller scattering matrix, 67, 69, 71 

Sj - scattering amplitude matrix, 68 

$ - target orientation angle, 22, 32 

$DF, 46, 74 

9 - target orientation angle, 22, 32 

Odf, 46, 74 

P - target orientation angle, 22, 32 
^DF, 46, 74 

ri - parameter for choosing scattering angles, 67 
?7 - parameter for selection of scattering angles, 66 
ai, 3.2 target vectors, 22 
(j)s - scattering angle, 61, 65 
9s - scattering angle, 61, 65 
tteff, 8, 21 

fij - amplitude scattering matrix, 68 
fij - scattering amplitude matrix, 68 
Qr - radiative torque efficiency vector, 30 
Qpr - radiative force efficiency vector, 30 

absorption efficiency factor (3abs> H 
ambient medium 

refractive index, 21 
averages over scattering angles, 64 

CALLTARGET, 47 
CALLTARGET and PBC, 48 
CMDFFT - specifying FFT method, 84 
CMDFRM - specifying scattering directions, 61, 84 
compiling and linking, 14 
conjugate gradient algorithm, 26 

ddscat.log_000 - output file, 25 

ddscat.out, 24 

ddscat.par 

ALLBIN, 25 

CSHAPE, 34 

diel.tab, 20 

DOTORQ, 30 

ETASCA, 30 

GAMMA, 21 

GKDLDR, 19, 28 

GPFAFT, 19 

lORTH, 67, 72 

IWRKSC, 25 

LATTDR, 19, 28, 29 

MXITER, 20 

NOTBIN, 19 



NRFLD, 20 
ORIBIN, 25 
PBCGST, 26 
PETRKP, 26 
SHPARi,SHPAR2,..., 34 
TOL, 26 

ddscat.par - BETAM1,BETAMX,NBETA, 33 
ddscat.par - PfflMIN,PHIMAX,NPHI, 33 
ddscat.par - THETMI,THETMX,NTHETA, 33 
ddscat.par parameter file, 18 
diel.tab - see ddscat.par, 20 
Dielectric Frame (DF), 37, 46, 74 
dielectric function of target material, 29 
dielectric medium, 1 1 
dielectric tensor, 73 

differential scattering cross section dCsca/dO,, 71 
DOTORQ, 30 

effective radius a^fi, 8, 21 

Electric field in or near the target, 75 

ELLIPSOID, 41 

error tolerance, 26 

ESELF, 81 

ETASCA, 30 

extinction efficiency factor Qext. H 

FFT algorithm, 27 
FFTW, 27 
Fortran compiler 

optimization, 15 

GAMMA - see ddscat.par, 21 
GKDLDR, 28 

GKDLDR - see ddscat.par, 19 
GPFA, 27 
GPFAFT, 54 

GPFAFT - see ddscat.par, 19 
GPFAPACK, 82 

IDVERR, 16 

IDVOUT, 16 

infinite cylinder, 50 

infinite hexagonal column, 54 

lORTH, 72 

lORTH parameter, 67 

iterative algorithm, 26 

IWRKSC - see ddscat.par, 25 

Lab Frame (LF), 31 
LATTDR, 28, 29 
LATTDR - see ddscat.par, 19 
LFRAME, 61 
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memory requirements, 30, 31 
Microsoft Windows, 13 
MKL 

Intel® MKL 
DFTI, 15 

MPI - Message Passing Interface, 16 

code execution, 24 
mtable - output file, 24, 25 

Mueller matrix for infinite targets, periodic in 1-d, 72 
Mueller matrix for infinite targets, periodic in 2-d, 72 
Mueller matrix for scattering, 67, 69, 71 
MXITER - see ddscat.par, 20 
MXNX,MXNY,MXNZ, 30 

NAMBIENT, 21 
NCOMP 

ddscat.par, 20 
nearfield calculation, 20 
NOTBIN - see ddscat.par, 19 
NRFLD, 75 

orientational averaging, 33 

nonrandomly-oriented targets, 34 

randomly-oriented targets, 34 
orientational sampling in /3, Q, and $,33 

PBC = periodic boundary conditions, 63 
PBCGST algorithm, 26 

periodic boundary conditions, 49, 50, 52-54, 56-59 

PETRKP algorithm, 26 

phase lag efficiency factor Qpha, 1 1 

PIM package, 26 

polarizabilities 

GKDLDR, 28 

LATTDR, 29 
polarization - elliptical, 68 
polarization of incident radiation, 64 
polarization of scattered radiation, 70, 71 
precision: single vs. double, 15 
PYD, 49, 53, 62, 63 
PZD, 49, 53, 63 

qtable - output file, 24, 25 
qtable2 - output file, 25 
qtable2 file, 25 

radiative force efficiency vector Qp^, 30 
radiative torque efficiency vector Qr, 30 
refractive index of target material, 29 
relative dielectric function, 1 1 
relative refractive index, 1 1 

scattering - angular averages, 64 
scattering angles 9s,(j)s, 65 
scattering by tilted cube, 64 



scattering directions, 61 
scattering efficiency factor Qsca. 1 1 
SCATTERLIB, 80 
size parameter x = kues, 8 
source code (downloading of), 13 
Stokes vector (/, Q, U, V), 69 

Target Frame, 73 
Target Frame (TF), 32 

target generation, 34 

target generation: infinite periodic targets, 48 
target orientation, 31, 32 
target routines: modifying, 47 
target shape options 

ANI_ELL_2, 38 

ANI_ELL_3, 39 

ANIELLIPS, 38 

ANIFTLPBC, 49 

ANIFRMFIL, 37 

ANIRCTNGL, 39 

BISLINPBC, 49 

CONELLIPS, 40 

CYLINDER1,40 

CYLNDRCAP, 40 

CYLNDRPBC, 49 

DSKBLYPBC, 50 

DSKRCTNGL, 40 

DSKRCTPBC, 52 

DW1996TAR, 41 

ELLIPSO_2, 42 

ELLIPS0_3, 43 

ELLIPSOID, 41 

FRMFTLPBC, 48 

FROM_FILE, 35 

HEX_PRISM, 43 

HEXGONPBC, 53 

LAYRDSLAB, 43 

LYRSLBPBC, 54 

MLTBLOCKS, 43 

RCTGL_PBC, 56 

RCTGLBLK3, 44 

RCTGLPRSM, 44 

RECRECPBC, 57 

SLAB_HOLE, 44 

SLBHOLPBC, 58 

SPH_ANI_N, 45 

SPHERES_N, 44 

SPHRN_PBC, 59 

SPHROID_2, 45 

TETRAHDRN, 47 

TRILYRPBC, 60 

TRNGLPRSM, 47 

UNIAXICYL, 47 
target.out - output file, 48 
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TFRAME, 61 
Vista, 13 

wOOOrOOOkOOO.pok - output file, 25 
wOOOrOOOkOOO.sca - output file, 25 
wOOOrOOOori.avg - output file, 25 
wOOOrOOOori.avg file - output file, 25 
Windows 7, 13 



